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A METHOD FOR THE THERMAL ANALYSIS OF 


SPACECRAFT, INCLUDING ALL MULTIPLE REFLECTIONS AND 
SHADING AMONG DIFFUSE, GRAY SURFACES 

By Doyle P. Swofford 
Langley Research Center 

SUMMARY 

A new method which uses finite surface elements has been developed and used to 
calculate temperature histories for spacecraft of arbitrary physical geometry. This 
method assumes gray surfaces with diffuse reflection and radiation properties and 
accounts for all multiple reflections. The analysis is performed in terms of the 
Cartesian coordinates of the four corners of each plane quadrilateral element. Shading, 
or optical blocking, is accounted for in the radiant heat exchange between each pair of 
surfaces and in the thermal fluxes from external sources which are incident on each 
surface. 

This method has been programed for a digital computer and the program can be 
used as is for the thermal analysis of any spacecraft whose surfaces may be approxi- 
mated by planar quadrilateral elements which obey Lambert 1 s cosine law. A listing of 
the program and its auxiliary programs is included in the report. 

An example problem of a complex, multiwinged earth-orbiting satellite is also 
presented. 


INTRODUCTION 

For a spacecraft which has extended members with large surface areas or has 
large enclosed spaces containing personnel or sensitive electronic equipment, shading 
and multiple reflections among its surfaces can be critical to its thermal design. 

The projected areas of the external surfaces of the spacecraft and the shape factors 
between all surfaces show the effect of shading of spacecraft surfaces by one another. 

The usual method of obtaining projected areas is by photographs of a model of the space- 
craft. Shape factors may also be obtained by a photographic method (ref. 1, pp. 399-402). 

For complex spacecraft in which temperatures are required at literally hundreds 
of nodal points, computer calculation methods are the only practical approach because of 
the large number of shape factors and projected areas that must be obtained. 



When the computer program in this report was written, other programs were 
already available which calculate projected areas and shape factors for configurations 
with shading and which account for multiple diffuse reflections of thermal radiation as 
well. (See ref. 2.) However, since these programs require much time to prepare and 
run, a more simple, yet flexible program, such as the one in the present investigation, 
was desired for development use in a wide variety of spacecraft. The thermal design 
may be fairly well fixed by use of the simpler program and then confirmed by use of a 
more sophisticated program. 

The treatment of the multiple reflections is made tractable by the usual approxi- 
mation that all the surfaces reflect and emit radiation diffusely (i.e., according to 
Lambert's cosine law). Once the shape factors and projected areas are calculated, the 
quadrature solution to this diffuse problem (developed by Gebhart, ref. 3) can be used. 
The basis of the approach used here is to divide the irregular quadrilateral surfaces of 
the spacecraft into grids by the use of two-dimensional coordinate systems embedded in 
each separate surface. The division allows the required integrations to be carried out 
in two-dimensional spaces. A shading test for each grid element is made by determining 
whether a line from the source point to the grid point is intercepted by any other surface 
belonging to the spacecraft. 

Calculated temperatures resulting from the shape factors and projected areas can 
be steady state or time dependent. If radiant heat transfer is predominant over conduc- 
tive heat transfer, equilibrium temperatures may be found by solving simultaneous equa- 
tions for the fourth power of the temperature of each surface. The nonequilibrium case 
(varying external and internal heat loads) is solved by integrating the time derivatives of 
the temperatures. In this case heat conduction is included in the calculation. 

SYMBOLS 


A area 

[A] diagonal matrix of the areas of a group of surfaces 

Aproj projected area of a surface 
[a] diagonal matrix of surface absorptivities 

B thermal flux incident on one surface from others 

{b} one-dimensional array of thermal flux incident on each surface from the 

other surfaces 
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specific heat of a material 


radiation coefficient between two identically shaped, closely spaced parallel 



matrix of the radiation coefficients Dy 
hemispherical emissivity of a surface 
diagonal matrix of the emissivities of a group of surfaces 
shape factor between surface i and surface j 

square matrix of shape factors between pairs of a group of surfaces 

shape factor of a spacecraft surface i for the albedo flux from a planet 
surface 

shape factor of a spacecraft surface i for planet- emitted radiation 
square matrix which yields { b \ when it operates upon {e} (see eq. (6)) 
sum of the thermal flux emitted by a surface and the flux reflected by it 
one -dimensional array of total thermal flux moving away from each surface 
conduction coefficient between nodes i and j 

square matrix of conduction coefficients between pairs of a group of nodes 
identity matrix 

shading indicator (1^ = 0: shaded by an intervening member; 

I s h = 1: not shaded) 

sunlight indicator (K = 1: sunlit; K = 0: shaded by planet) 


mass 



n 

q 

^cond 

<*int 

^rad 

r 


S 

T 



[W] 


O ,fty 


unit normal vector of a surface 

net heat flow rate to a node 

net rate of heat flow to a node by conduction 

rate of internal heat generation in a node 

net rate of heat flow between close parallel nodes by radiation 
position vector of a point 

position vector of point j relative to point i, r j - fj 
solar flux at position of planet or spacecraft 
absolute temperature, °K 

equivalent blackbody temperature of a planet, °K 
time 

position vector of the centroid of an elemental area 

position vector of the lower left corner of the jkth grid element of a plane 
quadrilateral 

unit vector directed toward the center of a planet from an orbiting spacecraft 
diagonal matrix with diagonal Wj = Dy 
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coefficient for heat conduction away from a node, 


z i = 


l 


h ij 


3 

diagonal matrix of coefficients for heat conduction away from each node 

parameters which each give the ratio of the lengths of two colinear vectors 
having a common origin 
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a',f3' abscissa and ordinate of a point in a normalized skewed coordinate system 

e sum of the thermal emission flux from a spacecraft surface and the reflected 

portion of incident thermal flux arriving at the surface directly from 
sources external to the spacecraft 

je } one -dimensional array of the fluxes e of the surfaces of a spacecraft 

0- angle formed by vectors nj and r jj 

9j angle formed by vectors nj and r 

9 S angle formed at the center of a planet by lines to the sun and to an orbiting 

spacecraft 

parameters giving any vector lying in a given plane as a linear combination 
of two given vectors in the plane 

M ratio of the projected area of a surface to its total area (unit projected area) 

{p} one -dimensional array of unit projected areas 

p reflectivity of a surface 

[p] diagonal matrix of reflectivities of surfaces 

Pp reflectivity, or albedo, of a planet surface for solar radiation 

a Stefan-Boltzmann constant 

total radiative heat flow from surface j which is incident directly upon 
surface i 

d$ji heat flux incident upon elemental area dAj which is received from dAj 

<£ a planet solar albedo flux incident upon a surface of a spacecraft 

$p planet thermal emission incident on a spacecraft surface 
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^rad 

4>s 


net thermal radiation flux received by a plate i from identical plates j 
parallel to and near it, Dj^T^ - 

solar radiation flux incident upon a spacecraft surface 




n 


CO 


function which is integrated to yield Fj^, 


( n i ' r l.l)( n i • r ti) 
”( f i) ■ f lj) 2 


half the angle subtended by a planet from the position of an orbiting 
spacecraft 


Subscripts: 
a albedo 

p planetary 

s solar 

t thermal radiation emitted by spacecraft surfaces or planet surface 

When an expression with an asterisk affixed as a superscript has a negative value, 
it is set equal to zero; positive values are unchanged. 

A bar over a symbol indicates a vector quantity. 

ANALYSIS 

Radiative Heat Balance 

The analysis of the heat transfer by diffuse emission and reflection of thermal radi- 
ation between one surface and the remainder of a group of surfaces is introduced by first 
considering just two surfaces i and j exchanging thermal radiation. The two surfaces 
are assumed to be isothermal. 

The heat flux incident upon the elemental area dAj of surface i which is 
received from dAj of surface j is given by 

Hi dAi 

dd > ii = — — ^ cos 9: cos 0. 

7T 2 1 J 

ij 
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and dAj, 

and 0j and 0j are the angles made with the normals to dAj and dAj, respectively, 
by the line r jj between them as shown in sketch A: 


where Hj is the heat flux leaving dAj, r^ is the distance between dAj 



Sketch A 


If r j is the position vector of dAj, and r j that of dAj, then 


cos 0j = n j • 



and 




where Hj and Hj are unit normal vectors to dAj and dAj. Thus, 


d *ji 




dAj 


Now, let 


(n t ■ ?,))(«, • f tj ) 







1J 
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so that 


d ^i = H 3^ dA 3 


If 


is the total heat flow incident upon surface 


i from surface j, then 


d 2 *jj = d*,j dA, - Hj^jj dAj OAj 


and 


Assuming Hj constant over surface j gives 

V H i <IA i " a j ‘ H /u A > A ; 

where the shape factor Fy is the average value of ^ over the two surfaces. Unlike 
the conventional shape factor, is not nondimensional and must be multiplied by Aj 

to yield the conventional quantity fjj. Also, Fj^ = Fy. The average flux incident on 
surface i from surface j is defined as 

Bi = = BijAjKj 

This is the flux incident on surface i from a single surface j. For more than one sur- 
face j, the total flux incident on i is given by the sum over j of the individual 
contributions: 

B rl f ii¥i (1) 

j 

If surface i is concave, Bj includes a contribution from surface i to itself. 
When Bj and are treated as the ith components of linear arrays |B} and {H}, 
it follows directly from equation (1) that 

{B}= [F][A]{H} (2) 

where [F] is a square matrix with elements Fjj, and [A] is a diagonal matrix whose 
diagonal Aj is the area of surface i. 
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The flux Hj leaving surface j is made up of ej (the emitted flux and the 
reflected portion of flux incident from sources external to the spacecraft) and the 
reflected part PjBj of the incident flux from other surfaces, where is the reflec- 
tivity of surface j. Thus, 


(3) 

and 6j includes the reflected part of any external radiation directly incident on j, such 
as sunlight. In the array form, 


= +[p]{B} (4) 

where [p] is a diagonal matrix of the reflectivities of the surfaces. 

Equations (2) and (4) are two matrix equations in the two arrays {B} and {H}; 
therefore, unique solutions may be obtained for them. Substituting equation (4) into 
equation (2) gives 


|b|=[f]W1 £ Mf][aM)b| 


(5) 


and, on solving for {B}, 


{Bf = [I - F Ap] _1 [FA]{e[ = [G]{e\ (6) 

where [I] is an identity matrix. Since [G] is a function only of the surface properties 
and geometry, it need only be evaluated once for each configuration. 

The net heat fluxes are given by the difference {b} - {h} between the incident flux 
and the flux leaving the surfaces. From equation (4), 

{b} - {H} = [I - p]{Bf -{e} 


or 


|B -H} = [a]{B} -M 


(7) 


( 8 ) 


When equation (7) is substituted into equation (6), the following equation is obtained: 

{B - H} = [aG - I]H 
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where [a] is a diagonal matrix of the absorptivities of the surfaces and the matrix [g] 
is given in equation (6). 

Equation (8) gives the net flux on each surface due to all emissions and reflections 
of thermal radiation from all the surfaces. Now, only the fluxes from external radiation 
sources which are directly incident on the surfaces remain to be accounted for. 

In the specific case of a spacecraft immersed in a real environment, there will be 
two regimes of thermal radiation flux - the principally short-wavelength solar flux and 
the principally long -wavelength emission from the spacecraft surfaces and from the 
planet surface. The external thermal radiation sources are direct and planet-reflected 
(albedo) solar flux and emission from the surface of a planet due to its temperature. 

Heat conduction as well as radiation will be accounted for. A diagram of the heat trans- 
fer for a single node is shown in sketch B: 


^s,i ^p,i 


°s,l 

B t,i 


Node i 




H s,i 

h m 


( q cond) i 


( q rad\ 


Sketch B 

If q^ is the net rate of heat input to member i, then the heat balance is given by 


5j_ 

Ai 


m i c i dT i 
A i dt 


( B s - H s)i * ( B t - H t)i - (*»**p), 


(“int), + Kond), 

Aj + (^radij 


( 9 ) 
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where 


m; 


mass of i 


specific heat of i 


absolute temperature of i 


<p„ i solar flux directly incident on i, both planet reflected and direct from the sun 

s,i 

cpp'i flux emitted from the planet surface and directly incident on i 

(^rad)- ne t ra diation flux to i from j, which is parallel to and near i 
(qint). rate of internal heat generation 
(q con d)- ne t ra * e conduction to i from other members 
Equation (8) gives for the solar spectrum 


- H sl - [s°8 - 


( 10 ) 


where e s ^ is the reflected part of the incident solar and albedo flux on i, that is, 
(p s 0 s )^. For long -wavelength flux 


}B t - H t }= [eG t - l]{etf (11) 

where is made up of reflected planet-emitted flux and the thermal emission of i. 

Thus 


e t,i = CTe i T i 


P tfi 4> 


P,i 


or in array form, 


hh CT [ e ]{' r4 } + [ I - e ]K} < 12 > 

where a is the Stefan-Boltzmann constant, [e] is a diagonal matrix of the emissivities 
of the surfaces, and {t^} is an array of the fourth powers of the temperatures. 


11 



Equation (9) becomes, in the array form, 

)B S - H s | + |B t - H t ) * |0 S + 0p) + + j^} + Kad} (13) 

Substituting equations (10), (11), and (12) into equation (13) gives 

= [ a sG s - I][Psl{<M + CT [ e G t - l][e]{T 4 }+ [eG t - l][l - e]{cfr p } + {<|> s + 

+ (¥j + M + K4 a*) 

and after simplification, 


{!}= Ps]P + GsPaltM + W[ ! + <¥" - e >]{*p} - "MP - G t e]{T 4 } 

'OB 

The net rate of heat conduction to member i is given by 

(“con*), - 1 h ii ( T i ' T i) = Z ( h ij T i) - ^ h iij T i 
where h« is the conduction coefficient between nodes i and j. In array form, 

i^condf = MW ~ MW = t h ' Z W ( 16 > 

where z is a diagonal matrix with z i = ^ hjj. 

i 

The net radiant heat to i from the identical panel j with which it forms an 
isolated enclosure is given by 


Kad^lMTj " T i 


where Dj^ is the radiation coefficient between j and i, given by 
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a 


D ii = T 


= D 


i _ 1 


ij 


Thus, 


i'Cradl = t D - "iM 


( 17 ) 


where 


W is a diagonal matrix with 


= Z D «- 


Substituting equations (16) and (17) into equation (15) gives the heat balance as 


{a} - + G s p s]{^s + 0af + C e ]\j- + G t^ ' e ^]{^p} 

- a[e][l - G t e]{T 4 }+ [D - W]{t 4 } 

+ [^][h - z]{T} + |^| (18) 


Only the external sources <p s and <p p ^ remain to be evaluated. The 


solar flux <£ s j 
the sunlit area of surface 


incident on i is the simplest, being proportional to the projection of 
i upon a plane perpendicular to a line to the sun (Aproj). : 


A i^s,i “ ( A proj) i KS 

< ^=(V) 1 KS=| v KS 

Or, in array form, 


j0 s [=KS{,i s [ 


(19) 


where K has the value 0 if the spacecraft is in the shadow of the planet and the value 1 
if not; and S is the local solar flux. 

The planet-emitted flux is given by 


<b . « f .oT 
r P,i P,i 


4 

P 
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and, in array form, 


{^p} ^pfp} 


( 20 ) 


where f„ , is the shape factor of surface i for the planet-emitted radiation, and T n 

P?- 1 - (A. ^ 

is the equivalent blackbody temperature of the planet krT p is equal to the average value 
over the planet of the emitted thermal radiation fluxY 


The planet solar albedo flux on surface i is given by 



or, in array form, 


l^al = P p S c° S *e s | £ a) < 21 > 

where p p is the mean reflectivity of the planet for sunlight, f ai is the shape factor 
of surface i for the planet albedo flux, 0 g is the angle formed at the center of the 
planet by the lines to the spacecraft and the sun, and cos*0 s = cos 9 S if cos 6 S is 
positive and equal to 0 if cos 9 S is negative. 

The factors f p and f a are difficult to calculate, even without the added compli- 
cation of partial shading of a surface. An approximation for them - exact when the panel 
is exposed to the entire portion of the planet surface which is visible from the position of 
the spacecraft - is 


fa,i * fp,i * («i • Vp)*sin 2 a; * #ip,i sin 2 co 

Here Vp is a unit vector pointing from the spacecraft toward the center of the 
planet, and co is the angle included by V p and a tangent from the spacecraft to the 
planet surface. The factor jUp^ approximates the effect of shading on f p> i and f a> i- 

When the relations for the external fluxes (eqs. (19) to (21)) are substituted into 
equation (18), the heat balance becomes 

KSM 1 + G s Ps ]{M + PpS cos*e s sin 2 w [a g ] [I + G s p s ]|Mp} 

+ oT p sin 2 co[e][l + G t (I - e)]|jUp| - <r[e][l - G t e]|T 4 } 

+ [D - W]{t 4 } + [i][h - z]{T} + j-^j (22) 
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For the equilibrium case, the net heat flow to each node is zero, and with negligible con- 
duction, equation (22) can be solved for {t 4 }: 

o{t 4 ^ = [e(l - G t e) + i(D - W)] |KS[a s ][I + G s p s ]|p s } + p p S cos*0 s sin 2 u>[a s ][I + G s p s ]|p p }. 

+ oT p sin 2 w[e] [i + G t (I - e)] |Pp} + j^J j (23) 

For a multifaceted node which is nearly isothermal, although its faces receive 
different thermal radiation fluxes, the total net heat flow rate is given by 


<li = 



where n is the number of faces of node i. For the nonequilibrium case, 


dT i_ Aj /q \ 
dt mjCjVA/j 


( 24 ) 


Equation (24) may be integrated numerically to yield the temperature history, beginning 
with given initial temperatures. 

Listings of the three computer programs which calculate shape factors, projected 
areas, and temperatures are shown in appendix A. They are written in the FORTRAN IV 
language. The shape factors and projected areas are run in separate programs because 
they need to be calculated only once for each configuration. In the main program, sur- 
face properties, materials, heat loads, attitude, and flight trajectory can be varied for a 
fixed configuration without recalculating the shape factors and projected areas. Punched- 
card outputs from the other two programs supply shape factor and projected-area inputs 
to the main program. 


Computation of Geometric Shape Factors and Projected Areas 

Routines were written to enable the digital computer to calculate [F] and {pj 
using as input the coordinates of the corners of each plane quadrilateral into which the 
spacecraft has been divided. The factor which makes the calculation of the projected 
area of a flat plate complex is the shading by intercepting surfaces. 

In order to calculate the unshaded projected area of a plane quadrilateral, it is 
divided into an n x n grid. Formulas for the coordinates of the grid points and areas 
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of the grid elements are developed in appendix B. The ratio of the projected area to the 
total area is given by 


^i " a- (^source ' *u) *sh,i dA i 


~ Aj ( V source ’ “i)^ I sh,i AA i 


( 25 ) 


where V gource is the unit vector in the direction of the source and 


*sh 


■C 


0 for shading 


for no shading 


The method of determining shading is developed in appendix C. 

In the computation of the shape factor between a pair of plane quadrilaterals, both 
are divided into n x n grids as before, and each element of one is paired successively 
with each element of the other. The shading test is performed by the same scheme as 
for the projected areas. The shape factor Fjj between surfaces i and j is given 
by the average (taken over both surfaces) of the function 


p « 4 U, " ij dAi ,IAj 

J i,j 

- 1 y -( fi i • f ij)( g i • f ij), 

AjAj L 4 


My AA i AA j 


1,1 


7Tr ij 


( 26 ) 


The position vector ry is given by r j - r where r^ and r j are the position 
vectors of the centroids of the elemental areas AA^ and AA^. The same formulas as 
for the projected-area routine are used for the area of AAj and its centroid r j. 

If a surface k is nonplanar and is approximated by a set of N planar surfaces, 
then is 
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( 27 ) 


* k = 



i=l 


N 

I 

i=l 


and Fjjj between the nonplanar surface and any other surface j is 

N 

£ AiF y 

Fkj = ^ ( 28 ) 

Z A i 

i=l 

In this investigation the nodal surfaces into which a spacecraft is divided are called 
panels. A large member may be divided into several pieces to conform more nearly to 
the assumption of constant temperature over each panel. Surfaces of the spacecraft which 
can intercept thermal radiation which otherwise would impinge on any of the panels are 
called shaders. There will be fewer shaders than panels if any of the plane surfaces is 
subdivided into more than one panel. 

EXAMPLE PROBLEM 

An example of the application of the computer program is the prediction of temper- 
atures on the proposed Meteoroid Technology Satellite. Figure 1 is a photograph of a 
model of the spacecraft shown attached to the last stage of its booster, with meteoroid- 
detector panels deployed. The cubical modules attached to the central structure are 
experiments for measuring velocities of meteoroids. The octagonal prism at the top is 
mainly solar cell area. In this example, the spacecraft is spinning about its axis of 
symmetry. 

The projected areas and shape factors are calculated in separate programs to be 
used later as inputs to the temperature -prediction program. The Cartesian coordinates 
of panel corners and shader corners are used as input for the projected-area and shape 
factor programs. Actual areas of the panels are also provided by the projected-area 
program as input to the main program. 

Direct inputs to the main program are absorptivities for solar radiation, hemi- 
spherical emissivities, the product of mass and specific heat for each panel, radiation 
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coefficients between back-to-back panels, and orbit and spacecraft-attitude parameters. 

A circular orbit with an altitude of 300 nautical miles (556 km) and an inclination of 38° 
to the equator was assumed. 

In this example, 104 nodes were used. No conduction is accounted for in this 
example. Because of the small thicknesses and large surface areas of the nodes, radia- 
tion will outweigh conduction by far. 

It was necessary to use a very small time increment (0.05 min) in the numerical 
integration of the temperatures. Some of the nodes consist of a sheet of plastic film 
0.00025 inch (0.00064 cm) thick, with a very small thermal capacity; therefore, the solu- 
tion is unstable for larger time increments. 

The temperatures all converged within two or three orbits, without oscillation, to a 
solution which repeated itself on subsequent orbits. The temperatures compare reason- 
ably with those obtained by a similar heat-transfer computer program which does not 
take multiple reflections into account. 

The temperature history is plotted for seven nodes, whose locations are indicated 
in figure 1. The temperature plot for the outer face of one of the upper velocity detectors 
(see fig. 2) follows the heat inputs to it very closely because it has very low thermal 
inertia. At the far left of the plot, the decline from maximum value of the earth albedo 
and earth thermal fluxes is noted, while the solar input remains constant. The computed 
value of the albedo and earth thermal flux declines to zero at about 24 min after perigee, 
after which the albedo remains zero and the earth thermal flux increases. At about 
31 min, the spacecraft enters the earth's shadow. During shadow, the earth thermal flux 
passes through a maximum and goes to zero again at about 71 min. Sunlight appears 
again at about 66 min. Albedo flux increases from zero at around 71 min to a maximum 
about midway through the sunlit period. The times of maximums and minimums in the 
thermal fluxes incident from external sources and in the temperature do not coincide 
exactly, since the node receives radiation from other parts of the spacecraft also. 

Figures 3 and 4 show the temperature variation over the orbit for a horizontal 
detector and a vertical detector, respectively. A schematic illustration of the cross 
section of these meteoroid detectors is given in figure 5, with the bumper shields on each 
side shown. 


CONCLUDING REMARKS 

A method has been developed for the thermal analysis of geometrically complicated 
spacecraft whose surfaces can be approximated by plane quadrilaterals with gray surfaces 
having diffuse reflection and radiation properties. Optical blocking between surface 
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elements is accounted for automatically. Multiple diffuse reflections are also accounted 
for. Listings of the computer programs which perform the calculations are included. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., May 1, 1970. 
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APPENDIX A 


COMPUTER PROGRAM LISTINGS 
Spacecraft Temperature Program 


PROGRAM ORBTEMP ( I NP.UT ♦OUTPUT ♦ TAPE5= I NPUT ♦ TAPE6=0UTPUT ) 

0090 

HEAT TRANSFER PROGRAM FOR DIFFUSE RADIATION ON MULT I PANELED SPACECPAFT080 

0070 

DIMENSION AR1 ( 104 ) , AR2 ( 104 ) *AP1 ( 104 ) ♦ APP1 (104) ,APP2 (104) ,CSH( 1 04 ) ♦ 

1ETAS ( 1 04 ) ♦ ET AE ( 1 04 ) ,FSP 1 (104,104) ,FSPP1 (104,104 ) ,FSPP2 ( 1 04 ♦ 1 04 ) ♦ 

2HCOND ( 104, 1 04 ) ♦ WATE ( 1 04 ) ,RRAD ( 1 04 ♦ 1 04 ) ♦ T ( 1 04 ) »XMU( 1 04 ) ♦ TET AS (20 ) ♦ 

3 YMU ( 1 04 ) ♦ TXMU ( 104,20) ,TTXMU(20 ) ,DT ( 1 04 ) ♦ THELAM (3 ,3 ) , URTH (3,3 ) , 

4 I PI VOT( 104), INDEX (104,2) 

DIMENS ION BTHETA(20 ) ,BAL PHA (20 ) , B5ET A (20 > , RM3PN (3,3 ) , RM 1 I (3,3), 

1RM30P (3,3) , RM3 A (3,3),RM1B(3,3) ,RM20 (3,3), RM3NTH (3,3) , RMOB ( 3,3 ) , 

2RMEB (3, 3 ) ,UVEC (3,104) , VEC ( 3 ) 

D I MENS I ON XMUSUN ( 42 , 7,19), TPH I S ( 36 ) , PH I S ( 36 ) , AX I S 1 ( 

1 3,50 ) , AXIS2 (3, BO ) ,TXMUS (360 ) ,VE1 (3 ) , VE2 (3 ) 

D I MENS I ON VE3 ( 3 ), RM I PH I (3,3), RMP I F ( 3 , 3 ) 

DIMENSION PHIE( 104 > 

14 FORMAT ( 1 2F6 * 2 ) 


15 FORMAT ( 5E 16*8) 0930 

16 FORMAT (6F1 3,3 > 0940 

17 FORMAT (2X,10E12#5) 0950 

18 FORMAT (//6H T I ME = E 1 5 * Q , 1 8X , 22H I NTERNAL TEMPERATURE =E15.6/> 0960 

1 9 FORMAT ( 1H1 , 9X , 9HT I ME ( M I N ) ♦ 7X , 1 4 HTHET A ( RAD I A NS ) ,5X, 1 5HALT I TUDE ( M I LE 097 0 

IS) ,5X, 15HH 0 =SHADE , 1 =SUN ♦ 9X , 3HAMP , 18X,7HETAS( 1 )//) 0980 

23 FORMAT! 1 HI , 1 OX, 14HVALUES OF UVEC// > 0990 

24 FORMAT ( 1 HI , 1 OX, 14HPR0GRAM INPUTS//) 1000 

25 FORMAT ( 2X , 6E20 • 6 ) 1010 


1984 FORMAT (3(2l4,E16»8) ) 

1 01 01 FORMAT ( 7F1 1 .8 ) 

1856 FORMAT (71 10) 

1905 FORMAT ( 8E 1 6 * 8 ) 

NT I ME= 1 

READ (5, 1856) NPANEL , NFAC , NHCOND ♦ NRRAD , N I NSEC , IFTLUP, I SP I N 
DO 3 1=1, NPANEL 

AR2 ( I )=0. 

APP2 ( I ) =0 # 

DO 3 J=l, NPANEL 
HCOND ( I ,J)=0# 

FSP 1 ( I , J )=0. 

FSPP1 ( 1 , J ) =0* 

3 FSPP2 ( I , J ) =0 • 

READ (5, 1 856 ) KEXT , K I NT , KRRAD , K HCOND 

READ (5,15 )RADE ,ALT, VELTO, VELR0,TI ,DTI ,TID,S 1030 

READ (5, 15)TE,BOLT,G,AU,RE 

1322 FORMAT( 1H058HSOLAR CONSTANT, S-B CONSTANT , EARTH TEMP* , EARTH REFLECT 
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1 I VITY) 

WR ITE<6* 1 322 ) 

WR I TE (6*17) S * BOLT « TE * RE 

READ (5* 15 ) ABET A t A I NCE * A I NCO * ALPHA * C OMEGA * DOME G A « OMEGAP * PH IN* THET AO 
IF< IFTLUP) 1 1 • 1 1 • 12 
12 READ ( 5 « 1 5 ) BALPHA 
READ ( 5 ' 1 5 ) BBET A 
READ (5* 15 )BTHETA 

11 READ ( 5 « 15) ( AR 1 ( J ) * J= 1 « NPANEL ) 

READ (5, 14) ( API ( J ) ♦ J= 1 .NPANEL) 

READ (5* 14) (APPl ( J ) * J = 1 * NPANEL ) 

READ ( 5 « 14 ) ( CSH ( J ) . J= 1 . NPANEL ) 

READ (5* 14 ) (WATEC J) • J=1 ♦ NPANEL ) 

WR ITE (6i 1313) 

1313 FORMAT < 1 HOI OHEXT* AREAS) 

WR 1 TE ( 6 t 1 7 ) AR 1 

1314 FORMAT ( 1H025HEXT* SURF. ABSORPT I V I T I ES ) 

WRITE (6* 1 314 ) 

WR I T E ( 6 • 1 7 ) API 

1316 FORMAT ( 1H023HEXT* SURF. EM I SS I V I T I ES ) 

WR I TE (6*1316) 

WR I TE ( 6 « 1 7 ) APPl 
1318 FORMAT ( 1 H07HWE I GHTS ) 

WRITE(6* 1 318 ) 

WR I TE (6*17) WATF 

1320 FORMAT ( 1 HO 1 4HSPEC 1 F I C HEATS) 

WR I TE ( 6 ♦ 1 320 ) 

WR I TE (6f 17) CSH 

READ (5. 14 ) CT( J ) ♦ J=1 . NPANEL ) 

IF (KEXT.EQ.O ) GO TO 1111 
DO 1492 K= 1 ♦ NF AC 

READ ( 5 . 1 984 ) I 1 * J 1 ,FSPP 1 ( I 1 t J1 ) t I2i J2 .FSPP1 ( I 2 ♦ J2 ) * I 3* J3 .FSPP1 ( I 3« 
1 J3 ) 

FSPP1 ( J 1 * I 1 ) =FSPP1 ( I 1 t J 1 ) 

FSPP1 ( J2i I 2 ) =F SPP 1 ( 12. J2) 

1492 FSPP1 ( J3» 13 )=FSPP1 ( 13. J3) 

DO 1776 1=1. NPANEL 

DO 1776 J=l. NPANEL 
FSPl ( I . J)=FSPP1 ( I • J)*AR1 ( J) 

FSPP1 ( I * J ) =FSP1 ( I • J ) 

RR AD ( I . J )=-FSPl ( I • J ) * ( 1 .-API ( J ) ) 

1 776 IF ( I .EO. J ) RRAD < I . J ) =RRAD ( I * J ) + 1 • 

CALL SIMEQ(RRAD. NPANEL. FSP1 « NPANEL . DETERM . IPI VOT . 1 04 . ISCALE) 

1861 FORMAT ( 1 H 1 60X6HFSP 1 -S/// ) 
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WR 1 TE (6.1 861 ) 

WR I TE (6*1 905 ) (( FSPl ( J .K ) .K = l .NPANEL ) . J=1 .NPANEL ) 

DO 1849 1=1. NPANEL 
DO 1849 J=l. NPANEL 

RRAD ( I , J)=-FSPP1 ( I, J)*( 1.-APP1 ( J) > 

1849 IF(I.EQ.J) RRAD ( I . J )=RRAD( I » J)+l . 

CALLS I MEQ (RRAD, NPANEL. FSPP1 .NPANEL .DETERM , I P I VOT . 1 04 , I SCALE) 

1862 FORMAT! 1 HI 60X7HFSPP1-S///) 

WRITE (6. 1862 ) 

WRITE (6 .1905) ( (FSPP1 (J,K).K=1 .NPANEL) ,J=1 .NPANEL) 

1111 CONTINUE 

RADO = RADE + AlT 
THETA = THETAO 
P = RADO*VELTO 
COSEGA=COS ( COMEGA ) 

S I NEGA=S I N ( COMEG A ) 

COlNCE=COS ( AINCE > 

S I NCE=S I N ( A I NCE ) 

IF (KINT.EO.O ) GO TO 1033 
READ ( 5 . 1 4 > ( AR2 ( J ) » J= 1 « NPANEL ) 

READ (5.14) (APP2< J) .J=l .NPANEL) 

DO 701 1 1=1 .NINSFC 

READ (5. 1984 ) I 1 . J 1 . FSPP2 ( I 1 . J 1 ) . 12. J2.FSPP2 ( 12. J2) . 13. J3.FSPP2( 13. 
1 J3 ) 

FSPP2 ( J 1 . I 1 ) =FSPP2 ( I 1 , J 1 ) 

FSPP2 ( J2 « 12 ) =FSPP2 ( 1 2 . J2 ) 

7011 FSPP2 ( J3 . 13 ) =FSPP2 ( I 3 . J3 ) 

DO 1929 1=1, NPANEL 
DO 1929 J=1 .NPANEL 
FSPP2 ( I . J)=FSPP2< I . J > *AR2 ( J ) 

RRAD ( 1 . J)=-FSPP2< I , J)* ( 1 •— APP2 ( J ) ) 

1929 IF(I.EQ.J) RRAD ( ! ♦ J > =RRAD (I « J ) + 1 • 

CALLSI MEG (RRAD, NPANEL. FSPP2, NPANEL. DETERM, IPI VOT. 104, ISCALE) 

1863 FORMAT ( 1 HI 60X7HFSPP2-S///) 

WRITE (6, 1863 ) 

WRITE (6. 1 905 ) ( (FSPP2 ( J .K ) , K=1 .NPANEL ) . J=1 .NPANEL) 

1033 IF(KHCOND.EQ.O) GO TO 1030 
DO 1010 1=1, NHCOND 

READ (5, 1984 ) I 1 , J1 , HCOND ( I 1 . J 1 ) . 12, J2.HCOND( 12. J2 ) . 13, J3,HCOND( 13. J 
13) 

HCOND ( J 1 , I 1 ) =HCOND ( 1 1 , J 1 ) 

HCOND (J2, 12 )=HCOND( 12. J2 ) 

HCOND (J3. 13) =HC0ND ( 13, J3) 

1010 CONTINUE 


1120 
1 130 
1 140 
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1324 FORMAT < 1 HI 23HCONDUCT I ON COEFFICIENTS) 

WRITE(6« 1 324 ) 

WR I TE I 6 ♦ I 7 ) HCOND 
1030 CONTINUE 

DO 207 I=ltNPANEL 
DO 207 J = 1 * NPANEL 
207 RRAD ( I ♦ J ) =0 # 

IF!KRRAD#EQ*0) GO TO 1020 
DO 1011 1 = 1* NRRAD 

READ ( 5 * 1984 ) I 1 * Jl *RRAD I I 1 ♦ J 1 ) ♦ I 2 ♦ J2 «RRAD ( 12 ♦ J2 ) ♦ I 3 ♦ J3 ♦ RRAD ( 13* J3 ) 
RRAD ( J1 t I 1 )=RRAD( I 1 ♦ Jl ) 

RRAD ( J2 t I 2 ) =RRAD ( 1 2 • J2 ) 

RRAD ( J3 ♦ I 3 ) = RRAD ( 13* J3) 

1011 CONTINUE 
1020 CONTINUE 

READ <5* 1856) NETASfNPHIS 
IF ( I SPIN) 1312«1312«1311 

1311 DO 420 I = 1 ♦ NPANEL 

420 READ (5* lOlOl ) (TXMU( NJ)iJ=l * NET AS) 

4891 FORMAT! 1 HOI 2HMU BAR TABLE) 

WR I TE ! 6 *489 1 ) 

WR I TE ! 6 ♦ 1 7 ) ( (TXMUt I ♦ J ) * J= 1 tNETAS ) *1 = 1 ♦ NPANEL > 

1312 READ (5# 14 ) (TETAS (M ) ,M=1 * NET AS ) 

4893 FORMAT! 1 HOI OHETAS TABLE) 

WRITE!6*4893 ) 

WR I TE (6*17) I TETAS! M) *M=1 * NET AS) 

DO 7 M=1 f NET AS 

7 TETAS !M ) = TET AS ! M )*• 0 1 74532925 

IF! ISPIN.EQ.1 ) GO TO 103 
532 DO 4 1=1* NPANEL 

DO 4 J = 1 * NET AS 

4 READ! 5# 101 01 ) IXMUSUN! I *J«K> ♦ !<=! *NPH1S) 


4895 

2849 

1701 

4897 

246 


FORMAT! 1 Hi 2QHMU! ETAS ♦PHIS) TABLES ) 

WR I TE ! 6 * 4895 ) 

DO 2849 1=1 ♦NPANEL 

DO 2849 J = 1 ♦ NET AS 

WRITE (6 ♦ 1 701 ) IXMUSUN! I,JiK)*K=l ♦ NPHIS > 
FORMAT !//2Xf 10El2#5) 

READ ( 5 ♦ 14) ITPHISINP) fNP=l ♦NPHIS) 
FORMAT! 1 HOI OHPH IS TABLE) 

WR I TE ! 6 ♦ 4897 ) 

WR I TE ( 6 ♦ 1 7 ) !TPHIS(NP) ,NP=1 ♦NPHIS) 

DO 246 M= 1 ♦ NPH I S 

TPH I S ! M ) =TPH I S ( M ) *• 0 1 745329252 
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531 READ <5* 16) < ( AX 1 S 1 < K ♦ J > *K= 1 * 3 > ♦ J= 1 ♦ NPANEL ) 

READ ( 5 ♦ 16) ( (AXIS2(K, J) ,K=1 *3) « J=1 ♦ NPANEL) 

READ (5* 16 ) ( CUVEC CK ♦ J ) * K» 1 * 3 ) « J = 1 « NPANEL ) 

DO 8 L=1 * NPANEL 
ST OR = AX I S 1 ( 1 *L ) 

AXIS1 (1 t L ) = AX I S 1 (2*L) 

AX I SI ( 2 ♦ L ) “AX I S 1 (3iL ) 

AXIS1 ( 3 ♦ L ) = STOR 
STOR = AX IS2 Cl tL ) 

AX IS2 ( 1 *L ) = AX I S2 ( 2 « L ) 

AXTS2 (2 *L ) =AXI S? (3, L ) 

AXIS2(3«L)=ST0R 
STOR=UVEC ( 1 *L) 

UVEC < 1 * L ) =UVEC ( 2 • L ) 

UVEC (2 »L )=UVEC ( 3« L ) 

6 UVEC ( 3 * L ) =ST OR 

C PROGRAM TO CALCULATE ENVIRONMENT OF SPACECRAFT 1610 

C 1620 

C 

C ROTMTX(M*N*ETA) COMPUTES THE ELEMENTS OF A ROTATION MATRIX M ABOUT AXIS N 
C (N=l TO 3) THROUGH ANGLE ETA 

C 

103 CALL ROTMTX ( RM3PN ♦ 3 * PH I N ) 

CALL ROTMTX ( RM1 I ♦ 1 • A I NCO ) 

CALL ROTMTX ( RM30P ♦ 3* OMEGA P ) 

CALL MULT ( RM I PH I * RM 1 I * RM3PN ) 

CALL MULT C RMP I F * RM30P iRMI PH I ) 

100 I F ( I FTLUP ) 102*102*101 

101 CALL FTLUPCTHETA* ABETA* 1 *20*BTHETA*BBETA ) 

CALL FTLUP ( THET A * ALPHA * 1 * 20 * BTHETA * BALPH A ) 

102 SINTHE = RADE/RADO 
COSTHE= ( 1 «-S INTHE**2 ) **#5 
OMEGA=T I *DOMFGA 
THET AN = -THET A 

CALL ROTMTX (RM3NTH* 3* THETAN ) 

CALL ROTMTX (RM3A* 3* ALPHA ) 

CALL ROTMTX ( RM 1 B * 1 * ABET A ) 

CALL ROTMTX (RM20* 2 • OMEGA ) 

MULT (MP* Ml *M2 ) MULTIPLIES 3X3 MATRICES Ml AND M2 * G I V I NG ( MP ) = ( M 1 ) ( M2 ) 

CALL MULT (RMEB*RM3A *RMP IF ) 

CALL MULT ( RMEB « RM 1 B « RMEB ) 2540 

CALL MULT ( RMEB * RM20 ♦ RMEB ) 2550 


1660 

1670 

1680 

1690 


24 



o r> 


APPENDIX A - Continued 


C RME3 IS COMPLETED 2560 

C ROTATION MATRIX RMEB ROTATES the INERTIAL REFERENCE FRAME (Z-AXIS coincident 
C WITH PLANET SPIN AXIS*X-AXIS PASSES THROGH THE VERNAL EQUINOX) TO COINCIDE 
c with the SPACECRAFT body axes 
C 

DO 170 J=li3 

1 70 VEC ( J ) =RMEB ( J ♦ 1 ) *CO SEGA + RMEB ( J ♦ 2 ) *S I NEGA*CO I NCE+RMEB ( J * 3 ) *S I NEGA* 

1 SINCE 

CALL MULT (RMOB*RM3A*RM3NTH ) 

CALL MULT ( RMOB * RM 1B» RMOB ) 

CALL MULT ( RM03 * RM20 ♦ RMOB ) 

ROTATION MATRIX RMOB GIVES ORBIT POSITION. COLUMN 1 GIVES COMPONENTS OF UNIT 
C VECTOR IN DIRECTION OF CRAFT FROM CENTER OF PLANET 
C 

IF ( ISPIN.EQ. 1 ) GO TO 1 

DO 160 1=1 *NPANEL 

DO 9 J= 1 * 3 

VE1 < J)=AXIS1 ( J* I ) 

9 VE2 (J )=AXIS2 ( J * I ) 

COSE=— UVEC ( 1 * I ) *RMOB ( 1 ♦ 1 > -UVEC (2*1) *RMOB ( 2 * 1 ) -UVEC (3»I) *RMOB (3*1) 

S I NE= ( 1 *-C0SE*C0SE)**.5 

I F ( COSE *LT * ♦ 1 E — 7 ) ETAE ( I ) = 1 • 57079633 

I F ( COSE • GE • • 1 E — 7 ) ETAE ( I ) = AT AN2 ( S I NE ♦ COSE ) 

IF ( (1 •— ABS(COSE) > • LT • 1 • E— 8 ) GO TO 1809 

COSE = — RMOB ( 1 * 1 )*VE1 ( 1 >-RMOB (2*1 )*VE 1 ( 2 ) -RMOB (3*1 )*VE1 (3 ) 

S I NE =— RMOB ( 1 *1 )*VE2( 1 ) -RMOB (2* 1 ) * VE2 ( 2 ) -RMOB ( 3 * 1 )*VE2<3> 

PH I E ( I )=ATAN2( SINE* COSE) 

GO TO 1810 

1809 PHIEC I ) =0 • 

1810 CONTINUE 

COSE=UVEC (1*1) * VEC ( 1 )+UVEC (2*1) *VEC (2 ) +UVEC (3*1 )*VEC ( 3 ) 

SINE* ( 1 •~COSE*COSE )***5 

I F ( COSE *LT • • 1 E— 7 ) ETAS ( I )= 1 . 57079633 

I F ( COSE .GE..1E-7) ETAS ( I ) =ATAN2 ( S I NE *COSE ) 

IF ( ( 1 *-ABS (COSE ) ) *LT. 1 .E-8 ) GO TO 1918 

COSE = VE 1 (1 ) * VEC ( 1 ) +VE 1 (2 )*VEC (2 )+VEl*(3 )*VEC (3 ) 

S I NE= VE2 ( 1 ) * VEC ( 1 ) +VE2 ( 2 ) *VEC ( 2 ) +VE2 ( 3 ) *VEC ( 3 ) 

PHIS ( I )=ATAN2( SINE* COSE) 

GO TO 160 
1918 PH I S ( I ) =0 • 

160 CONTINUE 2620 

GO TO 2 
1 CONTINUE 
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COSE=VEC (2 ) 

SINE = ( 1 ,~COSE*COSE)**,5 

IF( ABS(COSE) •LT i . IE-7 ) ETASUN=1 ,57079633 

IF ( ABS (COSE) ,GE, , IE-7 ) ETASUN= ATAN2 (S INE,COSE ) 

COSE=-RMOB (2,1) 

S I NE= ( 1 • -COSE*COSE > ** , 5 

IF(ABS(COSE),LT.* IE-7) ETARTH=1 ,57079633 

IF ( ABS (COSE ) • GE , • IE-7 ) ETARTH=ATAN2 (SI NE,COSE ) 

DO 6 1 = 1, NPA NEL 
ETAS ( I )=FTASUN 
6 ETAE( I )=ETARTH 
2 CONTINUE 

AMP=RM0B (1,1 )*VEC ( 1 )+RMOB (2 , 1 ) *VEC ( 2 )+RMOB (3, 1 )*VEC<3 ) 
H =1 ,0 

I F ( -AMP , GT , COS THE ) H = 0 , 0 
IF (AMP.LT.O.O) AMP=0,0 
C 


DO 1000 J=1,NPANEL 
I K = 0 

IF ( I SPIN) 1 357, 1 357, 1 358 

1357 DO 680 I=1,NETAS 
DO 680 K=I ,NPHIS 
IK= IK+1 

680 TXMUS ( IK)=XMUSUN( J, I ,K) 

NY=NETAS*NPHIS 

CALL D I SCOT ( ETAS ( J) ,PHIS( J) , TET AS , TXMUS , TPH IS, 22 , NY , NPH I S, XMU ( J ) ) 
CALL DISCOT (ETAF(J) , PH I E ( J ) , TET AS , TXMUS , TPH I S ♦ 22 , NY « NPH IS, YMU ( J ) ) 
GO TO 1000 

1358 DO 1001 1=1, NET AS 

1001 TTXMU ( I ) =TXMU ( J , I ) 

CALL FTLUP (ETAS ( J ) , XMU ( J ) * 2 , NET AS , TET AS , TTXMU ) 

CALL FTLUP (ETAE ( J ) , YMU ( J ) * 2 « NET AS , TET A S , TTXMU ) 

1000 CONTINUE 

* AT THIS POINT PROJECTED AREAS FOR SOLAR AND PARENT BODY ARE AVAILABLE 
DO 1002 I=1,NPANEL 

DT( I )=H*S*AR1 ( I )*XMU( I )*AP1 ( I ) +BOL T*TE**4*S I NTHE**2*AR 1 ( I ) *YMU ( I ) 
1*APP1 ( I > + S*RE*AMP*S INTHE**2*AR1 ( I )*YMU< I > *A P 1 ( I ) 


(DIR, SOLAR+DIR, Pt.ANET EMISSION* PLANET-REFLECTED SOLAR) 

C 

TEMP=0% 

DO 1 003 K = 1 ,NPANEL 

TEMP=TEMP+ ( APP 1 (K )*BOLT*T (K ) **4+B0LT*TE**4*S I NTHE**2* YMU <K ) * ( 1 ,- 
1 APP1 (K) ) ) *FSPP1 ( I ,K)*AR1 ( I ) * APP 1 ( I ) 


2730 
274 0 
2750 
2920 


2955 
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1003 CONTINUE 
C 

C (TEMP ADDS EXTERIOR LONG WAVELENGTH CONTR IB. TO PANEL I FROM ALL PANELS* 

C DIRECTLY AND BY ALL MULTIPLE DIFFUSE REFLECTIONS. CONSISTS OF 

C EXT. THERMAL EMISS. BY PANELS + REFLECTED PORTION OF PLANET EM 

C 

SFMP=0# 

DO 1004 K=1*NPANEL 

SEMP=SEMP+ (H*S*XMU ( K ) +S*RE*AMP*S I NTHE**2* YMU ( K ) )* ( 1 .-API (K ) ) 2975 

1 *FSP1 ( I *K >*AR1 ( I )*AP1 ( I ) 

1004 CONTINUE 

(SEMP ADDS SOLAR SPECTRUM CONTR IB. FROM ALL PANELS INCL. ALL REFLECTIONS 

CONSISTS OF PANEL— REFLECTED PORTION OF DIR. AND PLANET— REFL • SOLAR 

UEMP=0 . 

DO 1005 K=1,NPANEL 

UEMP=UEMP+B0LT*APP2 (K )*T <K ) **4 *FSPP2 ( I *K )*AR2 ( I )*APP2 ( I ) 2995 

1005 CONTINUE 

(UEMP ADDS INTERIOR LONG WAVELENGTH CONTR I B . TO I IF IT IS PART OF AN 

ENCLOSURE • CONSISTS OF EMISSIONS AND ALL REFL. FROM SURF. OF 
ENCLOSURE ) 

VEMP=0 • 

DO 1 006 K=1 .NPANEL 

VEMP=VEMP+HCOND (K,J)*(T(K)-T<I)) 

10C6 CONTINUE 

(VEMP IS THE CONDUCTION TO I ) 


1007 


WE MP = 0 % 

DO 1007 K=1*NPANEL 

WEMP±WEMP+RRAD ( K * I >* (T (K )**4-T < I ) **4 ) 
CONT I NUE 

WEMP=WEMP*AR1 ( I ) 


(WEMP GIVES RADIATION FROM ANOTHER PANEL PARALLEL TO I AND FORMING A FLAT 

ENCLOSURE WITH I. EXAMPLE IS FLAT PRESSURE CELL METEOROID DETEC. 


XEMP=-BOLT*T ( I ) **4* ( APP 1 ( I >*AR1 ( I )+APP2 ( I )*AR2 ( I ) ) 
CXEMP GIVES LOSS BY EMISSION FROM I ) 
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DT ( I ) =DT ( I )+TEMP+SEMP+UEMP+VEMP+WEMP+x£MP 
DT ( I ) =DT ( I )/WATE( I >/CSH ( I ) 

1002 CONTINUE 

DO 9843 1=1 t NPANEL 

9843 T( I )=T( I >+DT(I )*DTI 

I F ( MOD ( NT I ME ♦ 100 ) • EQ • 0 ) GO TO 113 

GO TO 114 3100 

1 1 3 CONTINUE 

WR I TE ( 6 « 1 7 ) ( T ( I > « I * 1 • NPANEL ) 

114 IF(TI«LT*TID) GO TO 99 3140 

GO TO 999 3150 

C 3160 

C ROUTINE TO CALCULATE RADIUS IN ELLIPTICAL ORBIT 3170 

C 3180 

99 CONTINUE 3190 

DDT I = • 1 *DT I 

DO 992 J = 1*10 3200 

ACCEL - (P**2)/RAD0**3 -G* ( RADE/R ADO ) **2 3220 

ACCEL 1 = ACCEL 3230 

DO 991 I = 1 »5 3250 

RADI =RAD0+VELR0*DDT 1+ ( 2 • *ACCEL + ACCEL 1 )*DDT I **2/6* 

VELR1 =VELR0+*5* (ACCEL + ACCEL 1 ) *DDT I 

991 ACCEL 1 = (P**2)/RAD1**3 - G* ( R ADF/RAD 1 ) **2 3290 

RADH = RAD1-RADE 3300 

ANVELO =P/RAD0**2 3310 

ANVEL1 = P/RAD 1 **2 3320 

THETA = THETA + ( ANVELO+ANVEL 1 > /Z • 0* ( DT I / 1 0 « 0 ) 3330 

RADO = RADI 3340 

992 VELRO = VELR1 3350 

RATIO * RADH/ALT 3360 

HIHT = RADH/5280 • 0 3380 

T I =DT I *FLOAT (NT I ME 5 

I F ( MOD (NT I ME *20 ) • NE « 0 ) GO TO 30 

WR I TE ( 6 * 25 ) T I « THETA ♦ H I HT * H * AMP * ET AS < 1 ) 3390 

30 NT IME=NTIME+1 

GO TO 100 3400 

999 STOP 34 10 

END 3420 

SUBROUTINE ROTMTX ( ROTOR * M * ROT ANG ) 

DIMENSION ROTOR (3*3) 

ROTOR ( M * M ) - 1 #0 
COSE = COS ( ROT ANG ) 

S I NE = S I N ( ROT ANG ) 

DO 1 I = 1 * 3 
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I F( I *EQ • M > GO TO 1 
ROTOR ( I ♦ M ) = * 0 
ROTOR ( M * I) = • 0 
ROTOR ( I * I )=COSE 

1 CONTINUE 

I F ( M * EQ • 1 ) GO TO 2 
I F ( M • EQ • 2 ) GO TO 3 
ROTOR ( 1 *2)=SlNE 
ROTOR (2*1 )=-SINE 
RETURN 

2 ROTOR ( 2 * 3 ) =S I NE 
ROTOR ( 3*2 ) =-SI NF 
RETURN 

3 ROTOR (1*3) = -S I NF 
ROT OR (3*1 )=SINE 
RETURN 

END 

3430 

SUBROUTINE ROTMTX CALCULATES THE ELEMENTS OF A ROTATION MATRIX*GIVFN THE 
AXIS NO* AND THE CCW ANGLE OF ROTATION 

SUBROUTINE MULT(A,B*C) 3450 

D I MENS ION A(3*3)*B(3,3)*C(3*3)*D(3*3) 

DO 1 1=1*3 3470 

DO 1 J=1 *3 3480 

D ( 1 , J)=0.0 

DO 1 K=1 *3 

1 D(I*J)=D(I*J>+B(I*K)*C<K*J) 

DO 2 1=1*3 


DO 2 J=1 *3 
2 A ( I * J)=D ( I * J ) 

RETURN 3550 

END 3560 


* 
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Definitions and Instructions for the Temperature Program 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

* 

* 

* 

* 

* 

c 

c 

c 

c 

c 

c 

* 

c 

* 

c 

c 

c 

c 

c 

c 

* 

* 

c 

c 

c 


c 

c 


heat transfer program for diffuse radiation on multipaneled SPA CEC RAF TO 80 

0090 

IN PROGRAM 0100 

ANGLE OUT OF ORBIT PLANE OF SPIN AXIS (YAW) 

ANGLE BETWEEN ECLIPTIC PLANE AND EQUATORIAL PLANE 0140 

ORBIT INCLINATION ANGLE 0160 

NEGATIVE ELEVATION ANGLE (PITCH) 0160 

ALTITUDE OF SPACECRAFT (FT) 0170 

FOR DAYLIGHT* AMP=COS ( ANGLE BETWEEN LINES TO THE SUN AND SPACECR 
AFT DRAWN FROM THE CENTER OF THE EARTH) FOR DARKNESS ♦ AMP =0 
ABSORPTIVITY OF EXTERNAL CAVITY TO P ENERGY , D I MENS 1 ONLESS 
ABSORPTIVITY OF EXTERNAL CAVITY TO PP ENERGY * D I MENS I ONLESS 
ABSORPTIVITY OF INTERNAL CAVITY TO PP ENERGY * D I MENS I ONLESS 
EXTERNAL CAVITY AREA, INCH SQUARE 
INTERNAL CAVITY AREA, INCH SQUARE 

DISTANCE FROM SUN (ASTRONOMICAL UNITS) 0200 

UNIT VECTOR IN PLANE OF PANEL 
CROSS PRODUCT OF UVEC AND AXlSl 

PITCH TABLE (20 PARTS) 0240 

YAW TABLE (20 PARTS) 0250 

ORBIT ANGLE TABLE (20 PARTS) 0260 

STEFAN-BOLTZMANN CONSTANT 

SUN ANGLE IN ELLIPTIC PLANE FROM VERNAL EQUINOX 0320 

SPECIFIC HEAT OF JTH COMPONENT , BTU/POUND/DEG 
TIME DERIVATIVE OF ORB I T ANGLE ( RAD I ANS/M I N ) 0340 

ORBIT ANGLE (RADIANS) 0350 

SPIN RATE OF SPACECRAFT ABOUT 2-AXIS 
TEMPERATURE INCREMENT OF THE ITH PANEL 0370 

TIME INCREMENT (MIN) 0380 

ANGLE OF SUN OUT OF ORBIT PLANE 0390 

PARENT BODY ASPECT ANGLE 
SOLAR ASPECT ANGLE 
THE SHAPE FACTORS INPUT TO THE PROGRAM ARE ONLY GEOMETRIC# THE FACTORS 
WHICH TAKE INTO ACCOUNT ALL REFLECTIONS AND ALL RE-REFLECTIONS ARE COMPUTED 
IN THE PROGRAM USING GEOMETRIC FACTORS AND THE REFLECTIVITIES OF EACH NODE# 
EXTERNAL SHAPE FACTOR ELEMENT FOR P ENERGY, /INCH SQUARE 
EXTERNAL SHAPE FACTOR ELEMENT FOR PP ENERGY, /INCH SQUARE 
INTERNAL SHAPE FACTOR ELEMENT FOR PP ENERGY, /INCH SQUARE 
ACCELRAT I ON O ■=■ GRAVITY (FT/MIN**2) 0420 

SUN I ND I C A TOR , D I MENS I ONLESS 
LINEAR HEAT CONDUCTION COEF , BTU/M I N/DEG 

INTFGER TO INDICATE USE OF BBETA , BTHETA , BALPHA 0490 

USE 0 FOR CONSTANTS AND 1 FOR TABLES 0500 

I SPIN IS ZERO FOR A STABILIZED VEH I CLE , NON- ZERO FOR SPINNING# 


SYMBOLS USED 
ABETA 
AINCE 
AINCO 
ALPHA 
ALT 

AMP 

API ( J) 

APP1 ( J ) 

APP2 ( J ) 

AR1 ( J) 

AR2 ( J ) 

AU 

AXIS1 

AXIS2 

BALPHA 

BBETA 

BTHETA 

BOLT 

COMEGA 
CSH ( J ) 

DDEL 

DEL 

DOMEGA 
DT ( I ) 
DTI 
EPI 
ETAE( J) 

ETAS ( J) 


* 

FSP1 (K, J ) 

* 

FSPP1 (K, J) 

■* 

FSPP2 ( K , J ) 

C 

G 

* 

H 

* 

HCOND ( K , J ) 

C 

IFTLUP 
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l 


< 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

* 

c 

c 

* 

c 

* 

* 

c 

c 

c 

c 

c 

c 

c 

c 

* 

* 

c 

* 

c 

c 

c 

c 

c 

c 


KEXT ♦ K 1 NT . KRRAD • AND KHCOND INDICATE WHETHER THERE ARE NON-ZERO ELEMENTS OF 
FSP1/FSPP1 ♦ FSPP2 ♦ RR AD . AND HCOND RESPECTIVELY 
NETAS NO* OF VALUES IN TABLE OF ETAS 

NF AC IS THE NO* OF CARDS USED TO LOAD THE NON-ZERO GEOMETRIC SHAPE FACTORS 
BETWEEN EXTERNAL SURFACES(3 TO A CARD)*FSPP1 IS INITIALLY USED AS A DUMMY 
NAME FOR THIS ARRAY* 

NHCOND IS SIMILARLY THE NO. OF CARDS FOR HCOND *NRRAD FOR RRAD .N I NSFC FOR THE 


INTERNAL GEOM. SHAPE FACTORS (READ IN AS FSPP2 > . 

NPANEL NUMBER OF PANELS 0510 

NPHIS NO. OF VALUES IN TABLE OF PHIS 

OMEGA ROLL ANGLE OF BODY 

OMEGAP ARGUMENT OF PERIGEE 0530 

P DENOTES SOLAR SPECTRUM . PP DENOTES THERMAL RADIATION SPECTRUM PEAKING A 
T MUCH LONGER WAVELENGTHS 

PH IN ARGUMENT OF ASCENDING NODE 0560 

PHIS(J) SOLAR AZIMUTH ANGLE 

RADE EARTH RADIUS (FT) 0590 

RADH ALTITUDE OF SPACECRAFT (FT) 0600 

RADO INJECTION RADIUS OF ELLIPTICAL ORBIT (FT) 0610 


REFE REFLECTIVITY OF PARENT BODY TO SOLAR ENERGY * D I MENS I ONLE SS 

RRAD(K.J) RADIATION COEFFF I C I ENT * W I TH NET HEAT TRANSFER RATE = 


AREA*RRAD (K. J >* (T (K )**4-T ( J >**4 ) 

S SOLAR CONSTANT. BTU/INCH SQUARE/M I NUTE 

SINTHE SINE OF THE EARTH VIEW ANGLE 0690 

T ( J ) TEMPERATURE OF JTH COMPONENT . DEG 

TE TEMPERATURE OF PARENT BODY. DEG 

TET AS . TPH IS. TXMU . TXMUS TABLES OF ET AS * PH I S . XMU 

THETA ORBIT ANGLE FROM PER I GEF 0720 

THETAO INITIAL THETA (ARBITRARY) 0730 

T I TIME (MIN) 0760 

TID CALCULATION TIME LIMIT 0770 

UVEC UNIT NORMAL VECTOR TO PANEL 

VELRO INJECTION VELOCITY NORMAL TO THE EARTH ( FT /M IN) 0810 

VELTO INJECTION VELOCITY TANGENTIAL TO THE E ARTH ( FT/M I N ) 0820 

WATE(J) MASS OF JTH COMPONENT . POUNDS 

XMU(J) UNIT PROJECTED AREA . D I MENS I ONLESS . SOLAR 


XMUSUN ( K ♦ I . J ) TABLE OF XMU OF EACH PANEL FOR GIVEN ETAS AND PH I S 

YMU(J) UNIT PROJECTED AREA ♦ D I MENS I ONLESS .PARENT BODY 

0850 


ARRAY DIMENSIONS FOR AN INDIVIDUAL CASE ARE- 

BTHETA.BALPHA.BBETA-SIZE OF TABLES OF PITCH AND YAW VERSUS ORBIT ANGLE. 
UVEC (3. NPANEL) 

AR1 .AR2.AP1 .APP1 . APP2 . CSH ♦ ETAS . ETAE- ALL ( NPANEL ) * 


* 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


FSP1 * FSPP 1 .FSPP 2 ♦ HCOND * RRAD * (NPANEL * NPANEL > • 

WATE * T ♦ XMU# YMU— ALL (NPANEL ) . 

TETAS ( NETAS ) ♦ TXMU (NPANEL ♦ NET AS ) tTTXMU ( NET AS ) *DT (NPANEL) . 

XMUSUN (NPANEL .NET AS ♦NPH IS ) ♦ TPH l S (NPH I S ) ♦PHIS (NPHIS ) • 

AXIS1 ( 3 ♦ NPANEL ) ♦ AX I S2 ( 3 ♦ NPANEL ) ♦ TXMUS (NET AS*NPH I S ) ♦ FAR ( NPANEL , NPANEL ) • 
I P I VOT ( NPANEL ) ♦ I NDEX ( NPANEL ♦ 2 ) . 


INPUT DATA LOADING OPDER 

1. (8110) NPANEL .NFAC ♦ NHCOND .NRRAD ♦ N I NSFC ♦ IFTLUP. I SPIN. 

2. (8110) KEXT.KlNT#KRRAD.KHCOND. 

3. (5E16.8) RADE ♦ ALT ♦ VELTO ♦ VELRO ♦ T I / DTI«T!D*S« 

4. (5E16.8 ) TE ♦ BOLT ♦ G ♦ AU«RE . 

5. (5E16.8 ) ABETA. AINCE *AlNCO. ALPHA ♦ COMEGA/DOMEGA ♦ OMEGAP ♦ PH I N , THET AO • 

6. IF ( IFTLUP.GT.O ) ♦ (5E16.8 ) B ALPHA ♦ BBET A ♦ BTHE TA • 

7. (5E16.8) AR1 8. (12F6.2) API 9 . ( 1 2F6 . 2 ) APP 1 1 0 • (1 2F6 • 2 ) CSH 

1 1 • ( 1 2F6.2 ) WATE 12.C12F6.2) T 

13. IF(KEXT.GT.O) FSPP1 (3 (2 14 , El 6.8 ) ) 

14. IF (KINT .GT.O ) (12F6.2) AR2/APP2 15. I F ( K I NT • GT • 0 > FSPP2 (3<2I4*E16.8 

16. IF(KHCOND.NF*0) HCOND ( 3 C 2 I A ♦ E 1 6 .8 ) ) 

17. I F ( KRRAD.NE . 0 ) RRAD ( 3 ( 2 I 4 ♦ E 1 6 . 8 ) ) 

18. (8110) NET AS ♦ NPH I S 

19. (7F11.8) TXMU 20. (12F6.2) TETAS 

IF ( ISPIN.EQ.O ) - 21. (7F11.8) XMUSUN^22. (12F6.2) TPHIS« 23. (6F13.3) 

AXIS1 ♦ AX I S2 ♦ UVEC 


PROGRAM OUTPUT - 

1. IF (KEXT.NE*0) ♦ (8E16.8 ) FSP 1 ♦ FSPP 1 -CONVERTED TO ACCOUNT FOR REFLECTIO 
2.IF (KINT.NE.O) ♦ (8E16.8) FSPP2-ALS0 CONVERTED 

3. EACH M I NUTE ♦ T I ( M I N. ) ♦ THET A ♦ ALT . (STAT.MI. ) ♦ H ( 1 • FOR SUNL I T ) ♦ A MP ( MEA - 

SURE OF PLANET-REFLECTED SUNL I GHT ) ♦ Et AS ( 1 ) • (2X.6E20.6) 

4. EVERY TEN M I NUTES ♦ TEMP • OF EACH PANEL ♦ DEG. R <2X*10E12.5) 


l 


i 
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Shape Factor Program 


PROGRAM SHAPE ( ! NPUT * OUTPUT * PUNCH • TAPE5= I NPUT * T APE6 = OUTPUT > 

C SHAPE FACTOR PROGRAM FOR DIFFUSE RADIATION 

C THE FOLLOWING INSTRUCTIONS MUST BE FOLLOWED WHEN USING THIS PROGRAM 
C IN THE MAIN PROGRAM CHANGE DIMENSION STATESMENT TO READ 

C DIMENSION C(M,4,3),D(N*5»3) •NUMB1 (M ) *NUMB2 (N) * COEFFA (L * L *4 ) t COEFFR (L«L*4 ) 

C IN SUBROUTINE _ ZAP CHANGE DIMENSION STATEMENT TO READ 
C DIMENSION Z ( N * 5 * 3 ) 

C IN SUBROUTINE AREA « CHANGE DIMENSIONS TO READ 
C DIMENSION COEFFA (L *L • 4 ) * COEFFR (L *L *4 ) 

C IN SUBROUTINE BLOCK CHANGE DIMENSION STATEMENT TO READ 
C DIMENSION B(N*5*3) *NUMB2(N) 

C WHERE M IS THE NUMBER OF PANELS AND N IS THE NUMBER OF BLOCKING PANELS.* L IS 
C THE GRID SIZE* 

C VARIABLE DEF I N I T I ONS * NO IS THE NUMBER OF PANELS • NBLOCK IS THE NUMBER 
C OF BLOCKING PANELS* C IS AN ARRAY CONTAINING THE COORDINATES OF PANELS* 

C D IS AN ARRAY CONTAINING THE COORDINATES OF THE BLOCKING PANELS 
C NUMB 1 IS AN ARRAY OF THE NUMBERS ASSIGNED TO THE PLANES OF EACH PANEL 
C NUMB 2 IS AN ARRAY OF THE NUMBERS ASSIGNED TO THE PLANES OF EACH BLOCKER 
C 

C PUNCHED OUTPUT - 

C (3(2I4*E!6*8) ) NON-ZERO SHAPE FACTORS FOR PANELS I AND J*3 TO A CARD 

C I 1 * J 1 *SHPFAC (IUJ1 )♦ 12* J2*SHPFAC (12* J2 ) * I 3 * J3*SHPFAC ( 13 • J3 ) 

C 

JD I MENS I ONCJ 1 00 * 4 • 3 ) ♦ D (75# 5* 3 ) * A ( 4 * 3 ) * B (_4 * 3 ) 

DIMENSION NUMB 1 ( 1 00 ) * NUMS2 ( 75 ) 

DIMENSION TENT (3) *NZ1 <3)_*NZ2<3> 

D I MENS I ON COEFFA ( 1 0 * 10*4 ) * COEFFR (10*10*4) 

COMMON/ONE/D 

■ COMMON /TWO /NUMB 2 
COMMON. /THREE/ _COEFFA_* CO.EFFR 

86 F ORMAT (7110) 

10 F0RMAT(3E16*8*2 !4 ) 

1 I FORMAT ( 1 2F6 • 1 ) 

330 FORMAT ’( 3 ( 2 I 4 * E 1 6 .8 ) ) 

READ ( 5 * 86 ) NO * NBLOCK 

. . _ READ <5.i 1 1 ) < ( (C t LvJ*K3 * K = 1 5 3.) *..J = l s.4JL *1 = 1 *.N0 ). 

READ (5* 11 ) t( (D( I*J«K)*K = 1 *3) *J=1 *4) *1 = 1 * NBLOCK ) 

READ (5*86 ) (NUMB 1 ( I > * I * 1 * NO ) 

READ (5* 86 ) (NUMB2 ( I ) * I = 1 * NBLOCK ) 

WR I TE (6*5) 

5 FORMAT ( 1 H 14HN0* OF P ANELS * 1 5HNO • OF BLOCKERS) 

_ WR I TE ( 6 * 86 ) NO* NBLOCK 
WR I TE ( 6 • T ) 

7 FORMAT ( 1 HOI 7HPANEL COORDINATES) 
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WRITE (6* I ! ) (((C(I#J#K)#K=1#3)#J=1#4)#I=1,N0) 

WR I TE (6*19) 

19 FORMAT ( 1 H 19HBL0CKER COORDINATES) 

WR I TE (6*11 ) {((D{I«JiK)»K=li3)iJ=h4)»I = l « NBLOCK ) 

WR I TE (6 # 6 ) 

6 FORMAT ( 1 HOI 6HPANEL PLANE NOS#) 

WRITE (6 #86) (NUMB 1 ( I ) # I = 1 # NO ) 

WR I TE ( 6 ♦ 8 ) * 

8 FORMAT ( 1 HOI 8HBL0CKER PLANE NOS#) 

WRITE (6 #86 ) (NUMB2C I ) # 1=1 #N5LOCK ) 

WRITE (6 #9) 

9 FORMAT (1 HI 4H I4H J3X9H ARF A OF I7X9HAREA OF J ) # 

NZ = 1 

NG = 5 

NG2=NG*NG 

F0URN2=FL0AT (4*NG2 ) 

TWON4= FLOAT ( 2*NG2*NG2 ) 

DO 75 M= 1 * NG 
I 1 =M-1 
121=11 +M 
N I =NG-M 
N I 1=2*N1+1 
DO 75 N = 1 * NG 
J1 =N-1 
J21=JM-N 
N J=NG— N 
NJ1=2*NJ+1 

COEFFA ( M «N # 1 ) = FLO AT ( N I *NJ+ ( N I + 1 )*(NJ+1 ) )/TW0N4 
COEFFA (M#N#2 )=FLOAT (II* (NJ+1 ) +M*N J ) /TW0N4 
COEFFA (M t N*3)=FL0AT( (NI+1 ) * J 1 +N I *N ) /TW0N4 
COEFFA(M#N#4)=FLOAT( I 1 * J 1 +M*N ) /TW0N4 
COEFFR ( M # N # 1 )=FLOAT(N! 1*NJ1 ) /FOURN2 
COEFFR ( M * N # 2 )=FLOAT( I21*NJ1 )/F0URN2 
COEFFR (M,N #3 )=FLOAT( I21*J21 )/F0URN2 
75 COEFFR ( M ♦ N * 4 ) = FLO AT (NT 1 * J2 1 )/FOURN2 
CALL ZAP (NBLOCK) 

IX=NO-I 
DO 3 1 = 1 * IX 
JX= 1+1 

DO 3 J=JX#NO 
DO 1 K=1 #3 
DO 1 L=l#4 
A (L#K)=C ( I «L ♦< ) 

1 B(L*K)=C(J#L#K) 
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NA = NUMB1 C I 1 

nb=numbi < J ) 

IF (NA.EQ.NB ) SAR=0*0 

IF (NA.FO.NB ) GO TO 4 

CALL FAKTOP ( A * B « $ AB ♦ SDA , SD8 * MA 4 MB # NG « MB LOCK 4 NG2 ) 

4 CONTINUE 

IFISAB.GT.l .F-10)G0 TO 333 

IF ( ( I .EG. IX) •AND* ( J*fq.N0> ) GO TO 333 

GO TO 3 

333 TENT (NZ ) = SAB 
NZ1 (NZ)=t 
NZ2 C NZ ) = J 

WP I TE ( 6 * 2 1 3 ) I 4 J*SAB 
213 FORMAT (2I4?E16*8*SH VIEWFAC) 

NZ “NZ+ 1 

1F(NZ*EQ*4> GO TO 33B 

I F ( ( I • EQ • I X ) • A NO • ( J • FQ • NO > ) GO TO 335 
GO TO 3 

335 PUNCH 330 * NZ 1 < 1 ) *NZ2( 1 ) ♦ TENT (1 ) *NZ1 (2 ) *NZ2 (2 ) * TENT (2 ) «NZ1 (3 ) *NZ2 (3 
1 ) 4 TENT < 3 ) 

NZ = 1 

3 CONTINUE 

STOP 
ENO 

SUBROUTINE ZAP(nBLOCK) 

D I MENS ION Z< 75 4543)40 ( 3) *C2 < 3 ) tC3 ( 3 ) ,C4 ( 3 ) ♦ C5 ( 3 ) ♦ AMI (3) 4 AM2 <3 ) 4 
1 V<3 ) 

COMMON/ONE/Z 
I =1 

5 DO 1 JO t 3 
C2(J)“Z(I <2 ♦ J)~7 ( I 4 1 i J) 

1 C3(J)=Z( I *4 4 J)-Z( I 4 1 4 J) 

CALL CR0SS<C5,C2#C3) 

CALL CROSS (Cl #C3*C5) 

CALL DOT(DET,Cl ,C2) 

CALL CROSS (AM2»C2iC5) 

DO 2 K = I t 3 
AMI (K ) =CI (K 5/DET 
2 AM2(<)=AM? (K ) /OFT 
DO 3 J=1 43 

3 V(J)=Z{M 4J) + 7(l4 3«J)-Z(I*24J)-Z(l44,J) 

C4(1)=AM 1 ( 1 )*V<1 )+AM 1(2 )*V ( 2 ) +AM 1(3)*V(3) 

C4(2)=AM 2(1 )*V(1)+AM 2C2)*V(2)+AM 2(3)*V(3) 

C4 (3 ) =0* 
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DO 6 J=1 #3 
Z( 1 *2. J)=C2 ( J) 

Z< I *3* J)*=C3( J) 

Z< I J)=C4< J) 

Z( r t5t j)=c5(j> 

1 = 1 + 1 

IF ( I »LT • ( NBLOCK+ 1 ) ) GO TO 5 

RETURN 

END 

SUBROUT I NE FAKTOR ( A ♦ B • S I J ♦ SDA * SDB • NA ♦ NB % NG , NBLOCK • NG2 ) 

DIMENSION A (4*3 ) «B (4«3 > *DA C 100 ) *DB ( 1 00 ) *PA ( 1 00*3 ) *PB ( 100*3 ) • AV<3 > • 
I BV < 3 ) • CV ( 3 ) • DV ( 3 > ♦ EV ( 3 ) • FV (3)*X(3)*Y(3) 

DIMENSION XA (3 ) ,CB (3) «RI J ( 16,3 ) * RDOTN ( 1 6 • 2 ) 4 GV ( 3 ) * HV ( 3 ) 

S I J = 0 • 0 
1 = 1 

DO 1 N= 1 #3 
A V ( N > = A (1 ,N) 

BV ( N ) = A (2 »N ) 

CV ( N ) = A ( 3 ♦ N ) 

DV ( N ) =B ( 1 »N) 

EV (N)=B (2«N ) 

FV ( N ) =B ( 3 ♦ N ) 

GV(N)=A(4*N) 

HV (N )=B (4 ,N ) 

R I J ( 1 *N)=DV(N)-AV(N ) 

RI J<2f N)=EV(N)-AV(N ) 

RI J(3«N)=FV(N)-AV(N) 

PI J (4«N)=HV (N)-AV(N) 

R I J ( 5 * N ) =DV ( N ) -B V ( N ) 

R I J ( 6 • N ) =EV ( N ) —R V ( N ) 

RIJ(7«N)=FV (N>— RVCN) 

R I J ( 8 • N ) =HV (N)-RV(N) 

RI J(9«N)=DV(N)-CV(N ) 

RIJ(10«N) =EV ( N ) — C V ( N ) 

R I J ( 1 1 * N ) =FV < N ) — CV ( N ) 

R I J ( 1 2 ♦ N ) =HV (N J — C V ( N ) 

R I J ( 1 3*N >=DV<N)-GV (N ) 

RIJ(14*N) =EV (N)-GV(N) 

RI J(15«N)=FV(N)-GV(N) 

RI J( 16*N)aHVCN)-GVCN) 

AV(N)=BV(N) — A V ( N ) 

BV(N)=CV(N)-BV(N) 

DV(N)=FV(N) — DV ( N ) 

1 EV(N)=FV(N)-FV(N) 
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CALL CROSS ( CV ♦ A V ♦ BV ) 

call cross (fv*dv*bv> 

DO 14 J=N16 
DO 1 4 N- 1 *2 

14 RDOTN ( J « N ) =0 • 

DO 1 5 J-l ♦ 16 
DO 13 N=1 *3 

RDOTN ( J * 1 ) =RDOTN C J ♦ 1 ) +R I J ( J ♦ N ) *CV ( N ) 

1 3 RDOTN ( J ♦ 2 ) =RDO TN ( J ♦ 2 ) +R I J ( J ♦ N ) *F V ( N ) 

I F ( ( RDOTN ( J # 1 ) • GE • 0 • 1 ) • AND • ( RDOTN ( J ♦ 2 ) • LE • — 0 • 1 ) ) GO TO 5 

15 CONTINUE 
SI J=C. 

SDA=1 . 

SDB = 1 • 

RETURN 

5 CALL DOT ( SA * CV ♦ CV ) 

CALL DOT(SB*FV»FV) 

SA=SA**.5 
SB=SB**#5 
DO 2 N= 1 « 3 
CV(N)=CV(N)/SA 

2 FV (N ) =FV(N ) /SB 

CALL AREA ( A *DA • RA t NG,NG2 ) 

CALL AREA ( B t DB • RR ♦ NG • NG2 ) 

DO 3 L=1 « NG2 
DO 3 J=1 * NG2 
DO 4 N=1 * 3 
X(N)=RA (L,N) 

X A (N)=X(N) 

Y<N)=RB( J,N) 

4 X(N)=X(N)-Y(N) 

CALL DOT (SI *CV*X> 

CALL DOT ( S2 « FV « X > 

I F ( ( S 1 • GE *0*0) • OR • ( S2 • LE* 0 • 0 ) ) GO TO 3 
CALL BLOKK (XA « Y «NBLOCK* G»NA , NB ♦ I) 

I F ( G • LT • • 5 > GO TO 3 
CALL DOT < S3 t X ♦ X ) 

S3=S3**2 

S I J=-S1 *S2/S3*0A ( L ) *DB ( J ) +S I J 

3 CONTINUE 
11 SDA=0 • 0 

SDB = 0 • 0 
DO 7 N=1 * NG2 
SDA=SDA+DA (N ) 
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7 SDB=SDB+DB<N> 

SI J=S! J/(SDA*SD«*3* 1415926536 ) 

IF ( ( ( S I J*SDA ) *LT* 1 5 ) • AND* ( ( S I J*SDB WLT • 1 *E-5 > ) SIJ=0« 

RETURN 

END 

SUBROUTINE AREA ( A ♦ DA • RA «NG « NG2 ) 

DIMENSION A ( 4 « 3 ) « DA ( I 00 )«RA(100*3) ♦ AV ( 3 ) ♦ BV (3)«CV(3) »DV ( 3 ) 

1 *T1 (3 ) *T2 <3 > *T3 (3 ) *T4 C3 ) tT5 (3) «T6 (3) «T7<3 > .T8 C3 ) 

D I MENS I ON COEFFA (10*10*4)* COEFFR (10*10*4) 

COMMON /THREE/ COEFFA * COEFFR 
18 FORMAT ( 4E 16*8) 

DO 2 N=1 « 3 
A V ( N ) =A ( 1 * N ) 

BV ( N ) = A (2«N) 

CV(N)=A (3*N) 

D V ( N ) = A ( 4 * N ) 

Tl (N)-BV(N>-AV(N) 

T2 (N)-CV(N) — BV ( N ) 

T3(N)=CV(N) — DV ( N ) 

2 T4 ( N ) =DV (N )-AV (M ) 

CALL CROSS ( T5 * T4 * T 1 ) 

CALL CROSS ( T6* T2 * T 1 ) 

CALL CROSS ( T7 * T4 * T3 ) 

CALL CROSS (T8*T2*T3) 

CALL DOT ( Z 1 ,T5* T5 ) 

CALL DOT ( Z2 * T6 ♦ T6 ) 

CALL D0T(Z3*T7,T7) 

CALL DOT ( Z4 « T8 * TS ) 

Z ! =Z1 ** • 5 
Z2-Z2**.5 
Z 3=Z3**.5 
Z4=Z4**.5 
K=1 

DO 1967 I ~ 1 « NG 
DO 1967 J=1,NG 

DA (K )=Z1 +COEFFA ( ! * J* 1 )+Z2*COEFFA ( I ♦ J ♦ 2 ) +Z3*C0EFFA { I * J ♦ 3 >+Z4*COEFFA 
l ( I • j*4 i 
DO 3 M= l * 3 

3 R A ( K ♦ M ) a A V ( M l^COE^FR ( ! , Jt 1 ) +BV ( M ) *COEFFR ( I , J • 2 ) + CV ( M ) *COEFFR < I * J *3 
l )+DV(M)*COEFFR< I ♦ J* 4 ) 

1967 K=K+1 
RETURN 
END 

SUBROUTINE BLOKK ( UV * VV * NBLOCK ♦ G * NU * NV * I ) 
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DIMENSION UV(3)*W(3>«BC75«5'3)tCI2<3)'CI3(3)*CI5(3) 

DIMENSION NUMB2 ( 75 ) « AM ( 3 ) * AM I 1 ( 3 ) , AM I 2 ( 3 ) « AM I 3 ( 3 5 * R I M I NA (3 ) 

COMMON/ONE/B 

COMMON/TWO/NUMB2 

EP I « I # OE — 06 

IJ = 1 

CALL SUB (AM* UV * VV ) 

3 CONTINUE 

IF ( (NU • EQ • NUMB 2 ( I ) ) • OR • (NV • EQ • NUMB 2 < I >) >GO TO 100 
DO 7 L=1 f 3 
CI2(L)=B( I • 2 «L ) 

C I 5 (L ) =B ( I *5*L) 

7 C! 3CL >*B< ! *3iL > 

CALL DOT ( DET *AM*CI5) 

IF (DET**2*LT.EPI) GO TO 100 
CALL CPOSS (AMU * A M « C 1 2 ) 

CALL CPOSS(AMI2 f CI3*AM) 

DO 8 M= 1 ♦ 3 
AMI1 ( M ) st AM I 1 (M) /DFT 
AMI2(M)=AMI2(M) /DET 
AM I 3 ( M ) =8 ( I ♦ 5 « M ) /DET 

8 R I M I NA ( M ) =UV ( M ) -B ( I * 1 « M > 

CALL DOT ( V3 ♦ AM I 3 ♦ R I MINA ) 

I F ( ( V 3 • GE *1 • 0 ) • OR • ( V 3 «LE.0*)> GO TO 100 
CALL DOTfVl * AMI 1 *RI MINA ) 

I F ( V 1 »LT *0*0) GO TO 100 
CALL DOT(Va*AMl?*RIMINA ) 

I F ( V 2 •LT.0*O) GO TO 100 

I F ( ( ( V 2 -1 #0)*(1 *0 + R( I «4«2> )-V 1 *B ( I * 4 * 1 ) ) • GT • 0 • 0 ) GO TO 100 
I F ( ( ( V 1 - 1 • 0 ) * ( 1 • 0+8 ( I « 4 • 1 ) ) — V 2 *B ( I *4 t 2 ) ) •GT«0«0 ) GO TO 100 
G = 0 • 0 
GO TO 60 
100 CONTINUE 

IF ( I J • EG * NBLOCK ) GO TO 98 
IF ( I • EQ • NBLOCK ) GO TO 10 
1 = 1 + 1 
GO TO I 1 
10 1 = 1 
11 IJsIJ+1 
GO TO 3 
98 G= 1 • 0 
60 CONTINUE 
RETURN 
END 
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SUBROUTINE SUB(C«A,B) 

DIMENSION C(3)«A(3)*B(3) 

DO 1 J=1 *3 
C( J)=A( J)-B(J) 

RETURN 

END 

SUBROUTINE CROSS ( C » A «B ) 
DIMENSION C (3 ) . A <3 ) »B (3 ) 

C ( 1 )=A (2 )*B <3)-R <2 )*A (3 ) 
C(2)=B(1 )*A(3>— A<1 )*B(3) 

C(3)=A( 1 )*B (2 )-B( 1 )*A (2 ) 

RETURN 

END 

SUBROUTINE DOT ( C ♦ A « B ) 

DIMENSION A (3) «R(3) 

C = A (I )*B ( 1 ) + A (2 >*B (2 )+A (3 >*8 (3 ) 

RETURN 

END 
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Projected Area Program 


PROGRAM PROJARC ! NPUT* OUTPUT * TAPES* I NPUT * TAPE6*OUTPUT ♦ PUNCH ) 

C PROGRAM FOR PROJECTED AREAS OF PLANE OU ADR I LATERALS * ACCOUNT I NG FOR SHADING 

C ARRAY DIMENSIONS WILL BE AS FOLLOWS- 

C MAIN PROGRAM- A ( NPNL * 4 * 3 ) • B (NSHDR , 5 « 3 ) « NUMP (NPNL ) * NUMS ( NSHDR ) * ETAS ( NM ) • 

C PHIS (NN) * RETAS (NM ) * RPH I S < NN ) « ARPAN (NPNL ) • XMUA V ( NPNL * NM ) * AND UNRMVEC ( NPNL * 3 ) • 
C SUBROUTINE PLANE— Z (NSHDR *5*3) 

C SUBROUTINE SH ADE-B ( NSHDR * 5 * 3 ) *NUMB2 CNSHDR ) 

C WHERE NPNL IS THE NUMBER OF PANELS * NSHDR IS THE NUMBER OF SHADERS *NM IS 

C THE NUMBER OF VALUES OF ETA—S* AND NN IS THE NUMBER OF VALUES OF PHI-S* 

C IN THE PROJECTED AREA PROGRAMING IS THE NUMBER OF EQUAL SEGMENTS INTO 

C WHICH EACH SIDE OF EACH PANEL IS DIVIDED TO FORM A GRID OF ELEMENTAL AREAS. 

C IF NG IS DIFFERENT FROM 1 0 • THE DIMENSIONS MUST BE DA ( NG**2 ) * RG (NG**2 , 3 ) « 

C COEFFA (NG*NG* 4 > * COEFFR ( NG ♦ NG ♦ 4 ) IN THE PROGRAM-AND IN THE SUBROUTINE AREA* 

C DA(NG**2 ) ,RA(NG**2«3> * COFFF A ( NG * NG * 4 ) .AND COEFFR ( NG * NG * 4 ) • 

C DEFINITIONS OF VAR I ABLES— 

C COEFFA (M*N,K ) IS ONE OF THE FOUR COEFFICIENTS IN THE FORMULA FOR THE AREA 
C OF GRID FLFMENT M * N OF A PANEL 

C COEFFR (M*N*K ) IS THE COEFFICIENT OF THE COORDINATES OF THE KTH CORNER OF THE 
C PANEL IN THE FORMULA FOR THE CENTROID OF ELEMENT M * N 

C NUMP IS AN ARRAY OF NUMBERS IDENTIFYING THE PLANE OE EACH PANEL 
C NUMS IS A SIMILAR ARRAY FOR THE SHADERS 

C A ( I < J*K )GI VES THE CARTESIAN COORDINATES OF CORNER J(J=1 TO 4> OF PANEL I 
C K=1 TO 3 CORRESPONDS TO X * Y * AND Z COORDINATES RESPECT I VELY« 

C B GIVES THE SHADER COORDINATES. 

C ETAS IS THE ARRAY ETA— S VALUES*PHIS OF THE PHt-S VALUES. 

DIMENSION A(50*4*3)*B(60*5*3) * NUMP (50) *NUMS (60) ♦ ETAS ( 37 ) 

1 * PH I S ( 72 ) * C ( 4 ♦ 3 ) * C 1 (3) *C2(3) *C3(3)* C4 ( 3 ) ♦ VNORM ( 3 > * V2 ( 3 > * SUN ( 3 ) * 

2RETAS ( 37 ) , RPH IS (72 ) • DA ( 1 00 ) ♦ RG ( 1 00 * 3 ) *RSUN(3) *RGR (3) 

D I MENS I ON COEFFA ( 10* ! O* 4 ) * COEFFR (10*10*4) 

DIMENSION ARPAN (50 )*XMUAV( 19 )♦ UNRMVEC ( 50 * 3 ) 

DIMENSION VEE 1(3) * VEF2 ( 3 ) 

D I MENS ION XMUPH I ( 72 ) 

DIMENSION AXIS1 ( 50 * 3 ) * AX I S2 (50 ♦ 3 ) 

DIMENSION C SET ( 37 ) *SNET(37> *SNPH(72) *CSPH(72) 

COMMON /ONE/B /TWO/NUMS 
COMMON /THREE/ COEFFA * COEFFR 

30 FORMAT (71 10) 

31 FORMAT ( 12F6. 1 ) 

READ (5*30 ) NPNL* NSHDR* I SPIN 
READ (5 *30) NG 

READ (5*30 ) (NUMP ( I ) * t * 1 * NPNL ) 

READ (5*30 ) (NUMS ( I > « 1 = 1 • NSHDR) 

READ (5 *31 ) ( ( (A( I « J*K) *K*1 *3) * J*1 *4 ) * 1 = 1 « NPNL ) 

READ (5*31 > ( ( (B( I * J.K) *K = 1 *3) • J=1 *4 ) * ! = 1 *NSHDR> 
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WR I TE ( 6 « 32 ) 

32 FORMAT (1H 10HNO. PANELS UHNO. SHADERS) 

WR I TE ( 6 ♦ 30 ) NPNL*NSHDR 
WR 1 TE ( 6 , 34 ) 

34 FORMAT ( 1 HO 1 2HS I ZE OF GRID) 

WR I TE (6 , 30 > NG 

WR I TE ( 6 « 35 ) 

35 FORMAT! 1 HOI 6HPANEL PLANE NOS.) 

WRITE (6. 30) (NUMP( I > . 1 = 1 «NPNL> 

WR ITE ( 6 « 36 ) 

36 FORMAT! 1 HOI 7HSHADER PLANE NOS.) 

WRITE (6.30) (NUMS! I ) . 1=1 .NSHDR) 

WR I TE ( 6 » 37 ) 

37 FORMAT! 1 HO 17HPANEL COORDINATES) 

97 FORMAT! 12F7. 1 ) 

WRITE (6 .97 ) ( ( ( A! I , J.K) ,K=1 .3) • J = 1 .4 ) . 1 = 1 ,NPNL ) 

WR ITE ( 6 * 38 ) 

38 FORMAT (1 HO 18HSHADER COORDINATES) 

WRITE (6 .97 ) ( ( (R ( I . J.K ) «K=1 ,3) . J=1 .4 ) , I =1 .NSHDR ) 
NG2=NG*NG 

FOURN2=FLOAT ( 4*NG2 ) 

TWON4 = FLOAT ( 2*NG2*NG2 > 

DO 75 M = 1 .NG 
I 1 =M-1 
I21=I1+M 
NI=NG-M 
N I 1 =2*N I + 1 
DO 75 N= 1 .NG 
J1 — N— 1 
J21 =J1 +N 
NJ=NG-N 
NJ1 =2*NJ+1 

COEFFA (M.N. 1 )=FLOAT ( N I *NJ+ ( N I + 1 )* (NJ+1 ) )/TW0N4 
COEFFA ( M.N. 2 )=FLOAT (II* (NJ+1 ) +M*N J ) /TWON4 
COEFFA (M.N .3 )=FLOAT ( (N I +1 ) * J 1 +N I *N ) /TW0N4 
COEFFA (M.N, 4 )=FLOAT ( I 1 * J 1 +M*N ) /TW0N4 
COEFFR ( M ,N . 1 )=FLOAT(NI 1*NJ1 J/FOURN2 
COEFFR(M,N,2 )=FLOAT( 121 *NJ1 ) /FOURN2 
COEFFR (M,N.3)=FLOAT( I21*J21 1/FOURN2 
75 COEFFR (M.N, 4 )=FLOAT !NI 1 *J21 )/FOURN2 
CALL PLANE (NSHDR) 

51 READ ( 5 » 30 ) NM.NN 
IF (EOF, 5 >998.999 
998 STOP 
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999 CONTINUE 

READ (5 t 31 ) (ETASf I) * 1 = 1 *NM ) 

READ ( *5 ♦ 31 ) (PH1S< I ) f Ial «NN ) 

WR I TE ( 6 « 33 ) 

33 FORMAT { 1 HOI 2HNO* OF ETA-S12HNO. OF PHI-S) 
WRITE (6 *30) NM « NN 
WRITE (6 *39) 

39 FORMAT ( 1 HOI 4HETA-SUN VALUES) 

WR I TE (6*31 ) (ET AS ( I ) * I = 1 » NM ) 

WRITE (6*40 ) 

40 FORMAT ( 1 HO I 4HPH ! — SUN VALUES) 

WRITE (6 *31 ) (PHI S ( I ) « 1=1 *NN) 

TOT=FLOAT(NN> 

ISHDR=1 

DO 150 LM=1 *NM 

RET AS (LM)=ETAS (LM)*1 . 745329252E-2 
CSET (LM ) =COS (RFTAS (LM ) ) 

150 SNET(LM)=SIN(RFTAS(LM) ) 

DO 1 75 MN=1 *NN 

RPHIS (MN)=PHIS (MN)*1 • 745329252E-2 
SNPH(MN)=SIN(RPHIS(MN) ) 

175 CSPH ( MN ) =COS CRPH I S < MN ) ) 

DO 1 IP=1,NPNL 
N I =NUMP ( IP ) 

DO 2 JA = I * 3 
DO 0 JC = 1 »4 

8 C(JC*JA)=A(IP« JC* JA) 

Cl ( JA )*C C2« JA )-C ( 1 • JA > 

C2(JA)=C(4* JA ) — C ( 1 « JA ) 

C3 ( JA ) 2=C (4 ♦ JA )— C(3*JA) 

2 C4(JA)=C(2*JA)-r(3,JA) 

CALL CROSS (VNORM, Cl *C2) 

CALL CROSS ( V2 « C3 « C4 ) 

CALL DOT < D * VNORM « VNORM ) 

CALL DOT ( D 1 *V2*V2) 

D=D**.5 
DO 68 JQ=1 « 3 

68 UNRM VEC ( I P * JQ ) = VNORM ( JQ ) /D 
D1=D1**%5 
ARMAC=#5* <D+D1 ) 

ARPAN ( IP) = ARMAC 

CALL AREA ( C , DA ♦ RG *NG , NG2 ) 

ARTOT = 0 % 

DO 13 J=M iNG2 
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13 ARTOT= ARTOT+DA ( J) 

WR I TE (6*16) IP 

16 FORMAT ( 1 HI 1 1 H PANEL NO.=!2) 

WR I TE (6*15) ARMAC 

15 FORMAT < 1 HO 2 OX ♦ SHARE A = F 1 6 • 0 ♦ 2X 9 7HSQ • IN* ) 

IF ( ( ABS (ARTOT/ARMAC-1 . ) ) *GE* *01 ) WR I TE ( 6 ♦ 1 4 ) 

14 FORMAT ( 1H037HARFAS DO NOT AGREE WITHIN ONE PERCENT) 

WR I TE ( 6 # 25 ) ARTOT 

25 FORMAT ( 1 H 23HSUM OF ELEMENTAL AREAS =El 6 * 892 X 9 7HS0* IN*) 

DO 3 KN=1 « 3 

3 VNORM (KNM VNORM (KN)/D 
DO 4 LM = 1 * NM 
SUM=0* 

WR I TE (6917) 

1 7 FORMAT ( t H03X4HETAS4X4HPH I S3X6HARPRO J I 0X2HMU ) 

DO 5 MN=1 *NN 

ARPRO J = 0 • 

IF ( ISPIN.EQ* 1 ) GO TO 27 
GO TO 21 

27 SUN ( 1 )»SNET(LM)*CSPH(MN) 

SUN ( 2 ) « SNET { LM ) *SNPH ( MN ) 

SUN ( 3 ) =CSET ( LM ) 

CALL D0T(SD0TN9VN0RM,SUN) 

IF (SOOTN.LE* 0. ) GO TO 23 
GO TO 47 
21 CONTINUE 

DO 20 K=1 *3 
29 VEE1 (< )=C1 (K ) 

SDOTN=CSET(LM) 

CALL DOTfVVl tVEFl *VEF1 ) 

VV1 =VV1 **.5 
DO 24 K=1 93 

24 VEE1 (K ) =VEE1 <K)/VV 1 

CALL CROSS ( VEE2 9 VNORM 9 VFE 1 ) 

DO 26 K=U3 

AXIS1 (IP 9 K ) = VEE 1 (<) 

AXIS 2 ( IP 9 K )=VEE2 (<> 

26 SUN ( < ) a VNORM (K ) *C SET ( LM ) +VEE MO *SNET (LM ) *CSPH ( MN )+VEE2 ( K ) *SNET ( LM 

1 ) *SNPH ( MN ) 

47 DO 7 N=1 *3 

7 SUN ( N ) = 1 00 0 • *SUN ( N ) 

DO 9 J=1 f NG2 
DO 10 K=1 93 

RSUN ( < )=RG ( J 9 K ) +SUN ( K ) 
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10 RGR (K ) = RG ( J , K ) 

CALL SHADE (RGR, RSUN 4 NSHDR , H 4 N I ♦ ISHDR ) 

I F ( H • LT • «5 ) GO TO 9 
ARPROJ=ARPRO J+DA ( J ) 

9 CONTINUE 

A RPRO J * A RPROJ * S DO TN 
23 XMU= ARPRO J/ARTOT 

WRITE (6# 1 8 ) ETAS < LM ) , PH I S ( MN ) « ARPROJ,XMU 
XMUPHI (MN) =XMU 
18 F0RMAT(2F8«1 ,2E16.8) 

IF(ISPIN) 5,04100 
1 00 SUM = SUM + A RPRO J 
5 CONTINUE 

TF ( I SPIN) 77»77,T8 

77 PUNCH 88, (XMUPHT (LP) ,LP = 1 , NN ) 

GO TO 4 

78 ARPJBR=SUM/TOT 
XMU8AR= ARP J0R/ARTOT 
XMUAV (LM ) =XMU0 AR 
WRITE(6,19) 

19 FORMAT ( 1 HOI 5HMFAN PROJ. ARE A • 7X * 1 OHA VFRAGE MU) 

WRITE (6 #20 ) ARP J3R * XMUB AR 

20 FORMAT ( 2E 1 6 • 8 ) 

4 CONTINUE 

IF( ISPIN#EG*0) GO TO 1 
PUNCH 88, (XMUAV(LM) ♦LM*1 ,NM ) 

1 CONTINUE 
88 FORMAT ( 7F1 1 .8) 

49 FORMAT ( 5E 1 6 • 9 ) 

PUNCH 49 , ( ARP AN (IP) t I P= 1 , NPNL ) 

50 FORMAT ( 6F 1 3 • 9 ) 

IF ( ISPIN) 72472,73 

72 PUNCH 50, ( (UNRMVEC ( IP, JQ ) 4 JO=l ,3 ) 4 IP= 1 ,NPNL ) 

PUNCH 50, ( (AXIS 1 C IP*K ) ,K=1 ,3) , IP-1 , NPNL ) 

PUNCH 50, ( ( AXIS2( IP*K) «K=1 , 3 ) , I P= 1 , NPNL ) 

73 PUNCH 30 4 NM,NN 

PUNCH 31 ♦ (ETAS(W) ,M = 1 , NM ) 

PUNCH 31 , (PHIS (N) ,N=1 ♦ NN ) 

GO TO 51 
END 

SUBROUTINE PL ANE ( NBLOCK ) 

DIMENSION Z(60*5.3)4C1 ( 3 ) 4 C2 ( 3 ) 4 C3 ( 3 ) , C4 ( 3 ) 4 C5 ( 3 ) 4 AM 1 (3) ,AM2 (3 ) 
1 V ( 3 ) 

COMMON/ONE/Z 
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1 = 1 

5 DO 1 J= 1 * 3 

C2C J)*ZC I «2« J)-Zt I* 1 * J) 

1 C3 ( J ) = Z ( I *4* J)-Z< I * 1 # J) 

CALL CROSS <CS*C2*C3) 

CALL CROSS (Cl * C3 * C5 ) 

CALL CROSS(AM2*C5«C2) 

CALL DOT { DET • C 1 *C2 ) 

DO 2 K=l»3 
AMI (K)=C1 (K > /DET 
2 AM2 (K > *AM2 (< )/DET 
DO 3 J=1 *3 

3 V(J)=2(MfJ) + Z(I«3*J)»Z<h2«J)-ZU»4«J) 

C4(l)sAM 1(1)*V(1)+AM 1(2) *V ( 2 ) +AM 1(3>*V<3) 

C4 (2 ) = AM 2(1 ) * V ( 1 ) + AM 2(2)*V(2)+AM 2(3>*V(3) 

C4 ( 3 ) =0 * 

DO 6 J = 1 <3 
Z( I * 2 * J ) =C2 ( J ) 

Z( I * 3 * J ) =C3 ( J ) 

Z( 1 * 4 * J ) =C4 ( J ) 

6 Z( I *5* J)=C5( J) 

1 = 1 + 1 

IF ( 1 »LT • (NBLOCK+l ) ) GO TO 5 
RETURN 
END 

SUBR • PLANS OPFRATES ON ARRAY Z (BODY COORDINATES OF THE 4 CORNERS OF EACH 
SHADER) TO FIND PARAMETERS NEEDED BY THE SHADING SURR» (SHADE) 

SUBROUTINE AREA ( A ♦ DA ♦ RA , NG % NG2 ) 

DIMENSION A <4*3 ) * DA ( 1 00 ) *RA ( 1 00*3) 4 AV<3 ) ♦ BV ( 3 > *CV(3 ) »DV(3) 

1 *T1 (3 ) * T2 (3 ) ,T3 (3 > • T4 (3 ) *T5 (3 ) *T6<3 ) * T7(3 ) ,T8 (3 ) 

D I MENS I ON COEFFA ( 10* 1 0*4 ) * COEFFR (10*10*4) 

COMMON /THREE/ COEFFA * COEFFR 

18 F ORM AT (4E16*8) 

DO 2 N= 1 * 3 
A V ( N ) = A ( 1 * N ) 

BV ( N ) = A ( 2 * N ) 

C V ( N ) = A (3,N) 

DV ( N ) = A (4*N) 

T 1 (N)=BV(N) — A V ( N ) 

T2 <N)=CV(N)-BV(N) 

T3(N)=CV(N)-DV(N) 

2 T4 (N)=DV(N)-AV(N) 

CALL CROSS ( T5 * T4 * T 1 ) 

CALL CROSS ( T6 # T2 * T 1 ) 
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CALL CROSS (T7.T4.T3) 

CALL CROSS (T8.T2.T3) 

CALL DOTtZl .T5.T5 ) 

CALL DOTIZ2.T6.T6) 

CALL D0T(Z3.T7,T7) 

CALL DOT(Z4,T8,T8) 

Zt=Zl**.S 

ze=z2**.s 

Z3=Z3**.<5 

Z4*Z4**.5 

K=1 

DO 1 967 I = 1 . NG 
DO 1967 J=1.NG 

DA (K )=Zl*COEFFA ( I « J. 1 >+Z2*COEFFA ( I . J . 2 ) +Z3*COEFFA ( I . J . 3 ) +Z4*C0EFFA 
1 ( I . J.4 ) 

DO 3 M=1 . 3 

3 RA(K,M)=AV(M)*COEFFR( I * J * 1 ) +8V ( M ) *COEFFR ( I , J « 2 )+CV ( M ) *COEFFR ( I .J.3 
1 ) +0 V ( M ) *COEFFR ( T , J . 4 ) 

1967 K=K+1 
RETURN 
END 

SUBR. AREA. USING THE COORDINATES OF EACH CORNER AND THE ARRAYS COEFFA AND 

COEFFR. COMPUTES THE COORDINATES OF THE CENTROID AND THE AREA OF EACH GRID 
ELEMENT OF A PANEL • 

SUBROUTINE SHADE ( UV * VV . NBLOCK , G . NU . I ) 

DIMENSION UV < 3 ) . VV (3) tB ( 60 .5.3) «CI2(3) * C I 3 ( 3 ) « C 1 5 ( 3 ) 

DIMENSION NUMB2 (60 ) « AM ( 3 ) • AM I 1 ( 3 ) « AM 1 2 ( 3 ) . AM I 3 (3 ) . AMI NX ( 3 ) 

COMMON/ONE/B 

COMMON/TWO/NUMB? 

EP 1=1 .0E-06 
I J = 1 

CALL SUB ( AM . VV . UV ) 

3 CONTINUE 

IF (NU.EO.NUMB2 ( I I )G0 TO 100 
DO 7 L— 1 . 3 
C I 2 <L > =8 ( I . 2 «L > 

C!5(L>=B(I.5.L> 

7 CI3(L)=B(I . 3 »L ) 

CALL DOT (DET » AM «C 15 ) 

IF(DET**2.LT.EP1 ) GO TO 100 
CALL CROSS (AMI 1 .AM. CI3 ) 

CALL CROSS (AMI2.CI2.AM) 

DO 8 M*1 .3 

AMU (M ) = AM I 1 (M ) /DET 
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AM 12 (M >=AMI2 <M)/DET 
AM 1 3 ( M ) =B ( I *5*M)/DET 
8 AM I NX ( M ) =B <1*1 ♦M)-UV(M) 

CALL DOT<V3*AMI3* AMINX) 

I F ( ( V 3 * GE • 1 • 0 ) • OR • ( V 3 #LE#0*>> GO TO 100 
CALL DOT (VI * AM I 1 ♦ AM I NX ) 

I F ( V 1 • LT *0*0) GO TO 100 

CALL DOT ( V2 * AM I 2 ♦ AM I NX ) 

I F ( V 2 #LT *0*0) GO TO 100 

IF ( ( (VI-1 « )* (1 #+B< I *4*2 >)-V2*B ( I *4* 1 ) )#GT#0. ) GO TO 100 
IF ( ( CV2-1 • >* <1 #+B ( I *4 * 1 > )-Vl*8 ( ! *4 *2 )) #GT#0# > GO TO 100 
G = 0 • 0 
GO TO 60 
100 CONTINUE 

IF ( I J#EQ#NBLOCK ) GO TO 98 
IF ( I .EQ.NBLOCK ) GO TO 10 
1 = 1 + 1 
GO TO 1 1 
10 1 = 1 
11 IJ=IJ+1 
GO TO 3 
98 G= 1 • 0 
60 CONTINUE 
RETURN 
END 

SUBR* SHADE DETERMINES WHETHER THE LINE BETWEEN TWO GIVEN POINTS INTERSECTS 
ANY OF THE SHADERS# 

SUBROUTINE SUB<C*4,B) 

DIMENSION A (3) *R(3> *C<3> 

DO 1 1=1*3 

1 C ( I ) = A ( I ) — B ( I ) 

RETURN 
END 

C SUBR# SUB GIVES VECTOR C=B-A 
SUBROUTINE CROSS ( C« A, B) 

DIMENSION A(3)*B(3) * C ( 3 ) 

C( 1 )=A(2)*B(3)«A(3)*B(2) 

C(2)=A(3)*B<1 ) —A ( 3 > *8 f 3 ) 

C(3)=A(1 )*B(2)-A(2)*B<1 ) 

RETURN 

END 

C SUBR# CROSS GIVES VECTOR C=CROSS PRODUCT OF VECTORS A AND B 
SUBROUTINE DOT(C«A#B) 

D r MENS ION A ( 3 ) « R ( 3 ) 


C = A ( 1 )*B<1 >+A(2)*B(2)+A (3>*B(3) 
rfturn 

END 

C SUBR# DOT GIVES C=DOT PRODUCT OF VECTORS A AND B# 
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DERIVATION OF CENTROID AND AREA FORMULAS 
FOR ELEMENTAL GRID SECTIONS OF A 
PLANE QUADRILATERAL 

The plane quadrilateral is divided into an n x n grid by dividing each side into 
n equal segments and connecting the corresponding dividing points on opposite sides. 
(See fig. 6 .) Counterclockwise from the lower left, the position vectors of the vertexes 
are A, B, C, and D. Any vertex may be taken as the starting point, so long as the 
order is counterclockwise. This gives the proper sense to the computed normal vector 
to the plane. 

The lower left corner of each elemental area is denoted by (j,k), where j and k 
are integers increasing from 1 to n; j, from left to right; and k, from bottom to top. 
Now, each segment of the line between A and B represents a change in position 
B A 

of — - — . Thus, the position vector of the point (j,l) is given by 
n 

V jfl = A+ (j - l)B-I-A=I[(n -3 + 1)A+ (j - 1 )b] (Bl) 

Similarly, along the line from D to C (where k = n + 1), 

V + 1 ' ff[ (n - 3 + 1)6 + <J - 1)5] (B2) 

The position vectors of the points (l,k) and (n+l,k) correspondingly will be given by 

v l k = i[( n - k + 1)A + (k - 1)d] (B3) 

Vi,k - M (n - k * i > 6 *( k - « 5 ] < b4 > 

If a line is drawn between Vj ^ and Vj n+ ^ and another is drawn between ^ and 

V n+ i k , their intersection will be Vj k , the position vector of the lower left corner of 

element jk. The line between Vj and Vj n+ j is divided by the intersection point 

in the same ratio as sides BC and AD are divided. Thus, in terms of V- •, and 
__ 1 > 

Vj n+ 1 ’ V j k may be expressed as 
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APPENDIX B - Continued 


V = s[("- k+1 >V +<k - 1)f i,n + l] 

Substituting equations (Bl) and (B2) into equation (B5) yields 

V- k = -~[(n - j + l)(n - k + 1)A + (j - l)(n - k + 1)B 
+ (j - l)(k - 1)C + (n - j + l)(k - 1 )d] 


(B5) 


(B6) 


This is the exact value of the position vector of the lower left corner of elemental area 
jk of the quadrilateral. The position vector of the centroid of the elemental quadrilateral 
may, except for extremely irregular shapes, be well approximated by 



This is the mean of the position vectors of the four corners of the grid element and would 
also be the point of intersection if the grid were twice as fine. 

The area of the element jk is found by dividing it into two triangles and taking 
half the magnitude of the cross product of two sides of each triangle. Since the upper 
right triangle of element jk is congruent to the lower left triangle of element j+1, k+1, 
the area of element jk is given by 


(AA) j,k " 2 


V j+l,k - V j,k) x ( V j,k+1 " V j,k) 


+ ! 


V j+2,k+l ' V j+l,k+l) 


X ^ V j+l,k+2 V j+l,k+lj 


(B8) 


Substituting equation (B6) into equation (B8) yields 
(AA)j,k = |[(n - k + 1)(B - A) + (k - 1)(C - D)] x [(n - j + 1)(D - A) + (j - 1)(C - B)] 

+ [(n - k)(B - A) + k(C - D)j x [(n - j)(D - A) + j(C - B)]| 


(B9) 
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APPENDIX B - Concluded 


Expanding equation (B9) gives 

(AA). k = j[(n - j + l)(n - k + 1) + (n - j)(n - k)j |(B - A) X (D - A) 
+ [(j - l)(n - k + 1) + j(n - k)] |(B - A) x (C - B)| 

+ [(n - j + l)(k - 1) + (n - j)k]|(C - D) x (D - A)| 

+ [(j - l)(k - 1) + jk]|(C - D) x (C - B)|| 


(BIO) 
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APPENDIX C 


DETERMINATION OF SHADING 

Five tests are made in the computer program to determine whether points X 
and Y are shaded from one another by the planar quadrilateral ABCD: 

C 


B 


The first test determines whether the point of intersection P of the plane of ABCD by 
the line XY falls between points X and Y. Then for each of the four sides of ABCD, 
it is found whether P lies toward the inside or the outside of ABCD from that side. 

The position vector of a point on the plane of the quadrilateral ABCD may be 
expressed as a linear combination of two vectors in the plane added to the position vector 
of a corner, say corner A: 



r plane = A + a'(B - A) + /3’(D - A) 

If P is the position vector of the point of intersection of the plane by the line con- 
necting X and Y, then 

P - X = y(Y - X) 

Setting P = ? plane gives 

a'( A - B) + /3’(A - £>) + y(Y - X) = (A - X) 
or 

a'V 1 + /3’V 2 + yV 3 = V 4 

This vector equation is, of course, three simultaneous linear equations in a ’, /3', 

and y - one for each vector component. 
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APPENDIX C — Continued 


The solution for y will be carried out first, since the quadrilateral ABCD will 
be immediately eliminated as a shader if the condition 0 < y < 1 is not met: 


( V 0 X 

( V 2) x 

( y 4) x 


( v l) Y 

( V 2) y 

( V 4) y 


( v l) z 

( V 2) z 

( V 4) z 

_ y i ' ( y 2 X y 4 ) 

( v i)x 

( V 2) x 

( V 3) x 

y l ' ( y 2 X V 3 ) 

( v l) Y 

( V 2) y 

( V 3) y 


( v *) z 

( V 2) z 

( y 3) z 



If a' < 0, the point lies to the left of side AD and if /3' < 0, the point lies below 
side AB. These will also be evaluated one at a time: 

a ,_V]V^3) Vfaxvs) 

Vj ■ (V 2 X V 3 ) V, • (V 2 X V 3 ) 

If a' and (3' are both greater than zero, more testing is required. 

The reason for the primes is that a' and j8' are actually quadratics in the two 
linear parameters for the skewed coordinate frame formed by the quadrilateral. Let a 
be a linear parameter characterizing a point moving from A to B or D to C 
as a varies from 0 to 1. Let j3 be the parameter for points along AD or BC as 
shown in the following sketch: 
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APPENDIX C - Continued 


Then, 

P = A + a(B - A) + /g|[D + a(C - D)] - [X + ce(B - A)]} 

= A + a(B - A) + /3(D - A) + a\ 3(A-B+C-D) 

The vector (A-B+C-D) lies in the plane of the quadrilateral and thus can be given as a 
linear combination of (B - A) and (D - X): (A-B+C-D) = X^(B - A) + X 2 (D - X). The 
parameters Xj and X 2 may be found by use of the vectors reciprocal to (B - X) and 
(D - A), denoted by superscript R and characterized by 

(B - X) R • (B - X) = 1 
(B - X) R • (D - X) = 0 
(D - X) R • (B - X) = 0 
(D - X) R • (D - X) = 1 
where, if N = (B - A) x (D - A), then 


( B-A) R =r 

I (D - A) X N I • (B - A) 

(S5- Rx (5- A) 

N X (B - S) • (D - A) 


Now, 


(B - A) R • (A-B+C-D) = X 1 


and 


(D - A) R • (A-B+C-D) = X 2 


Thus, 


P = A + a(B - A) + /3(D - A) + o^jx^B - A) + X 2 (D - A)j 
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APPENDIX C - Continued 


or 


P = A + a^l + Aj/^B - A) + 0^1 + AgCe^D - A) 


and 


a' = (1 + A^a 

p = (i + v *)' 3 

from which the values of a and 0 can be found directly. 

An advantageous alternative to finding a and 0 is found by plotting the quadri- 
lateral ABCD in a' 3' coordinates: 


Q {<*',?) 


C(1 + Aj, 1 + Ag j 



P(a',P) 


The a’ and /3' coordinates of the corners are found from the two equations for a’ 
and /3' in terms of a and 0, with a and 0 having a value of 0 or 1 at the corners. 

The points P and Q represent possible points of intersection between the plane 
of ABCD and a line connecting two points of interest. If point P lies on the right- 
hand side of line BC, then the single nonzero component of the cross product of 
vector (P - B) upon (C - B) will be positive. Thus for no shading, this third component 
will be given by 


[(P - B) X (C - B)] 3 = ( a ' - 1)(1 + A 2 ) - 0'Aj > 0 

If this quantity is zero or negative, P lies on or to the left of line BC. Similarly, the 
condition that the point lie above line CD is 

[(C - D) x (Q - D)] 3 = (1 + Aj) (0* - 1) - A 2 o-' > o 
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APPENDIX C - Concluded 


The five tests for shading, in the order of execution, are thus 

(1) 0 < y < 1 

( 2 ) a' % 0 

(3) ^0 

(4) (, a ' - 1)(1 + X 2 ) - j3’ • X x SO 

(5) O' - 1) (1 + X x ) - a' • X 2 = 0 

If any one of these fails, there is no shading between the two points by the quadrilateral 
ABCD. They must all hold for shading to occur. 
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Figure 5.- Schematic cross section of a pressurized-cell meteoroid detector with bumpers. 
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Figure 6.- Plane quadrilateral with n x n grid formed by dividing each side into n equal segments, with coordinate parameters j 

and k illustrated. 
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